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Abstract 

Computations are presented for one-dimensional, 
strong shock waves that are typical of those that form 
in front of a reentering spacecraft. The fluid mechan- 
ics and thermochemistry are modeled using two dif- 
ferent approaches. The first employs traditional con- 
tinuum techniques in solving the Navier-Stokes equa- 
tions. The second approach employs a particle simu- 
lation technique (the direct simulation Monte Carlo 
method, DSMC). The thermo chemical models em- 
ployed in these two techniques are quite different. The 
present investigation presents an evaluation of thcr- 
mochemical models for nitrogen under hypersonic flow 
conditions. Four separate cases are considered that are 
dominated in turn by vibrational relaxation, weak dis- 
sociation, strong dissociation and weak ionization. In 
near-continuum, hypersonic flow, the nonequilibrium 
thermo chemical models employed in continuum and 
particle simulations produce nearly identical solutions. 
Further, the two approaches are evaluated success- 
fully against available experimental data for weakly 
and strongly dissociating flows. 


Introduction 

A space-vehicle passing through the earth*s atmo- 
sphere will traverse a number of different flow regimes. 
At lower altitudes, the fluid density is sufficiently large 
for the flow to be considered in thermo chemical equi- 
librium. However, as the vehicle ascends higher into 
the atmosphere, the molecular collision rate will fall, 
and low-density effects will become increasingly im- 
portant. 

Continuum methods are successfully applied to 
flows in which the collision rate of the gas is suffi- 
cient to maintain Boltzmann energy distributions for 
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the various thermal modes of the gas. It is not neces- 
sary that the temperatures associated with each of the 
different modes be equal, or that chemical equilibrium 
prevails. Particle methods, such as the direct simu- 
lation Monte Carlo method (DSMC), are successfully 
applied to flows in which a reduced collision rate no 
longer supports equilibrium energy distributions. As 
the numerical cost of this technique is proportional to 
the fluid density, application has been limited to rar- 
efied flows. 

The computation of flow properties for the flight 
trajectories of many space vehicles will require the 
use of both continuum and particle methods men- 
tioned above. The interface between the different flow 
regimes is therefore of great importance. Clearly, it is 
desirable to obtain consistent results with these numer- 
ical methods in an overlapping near-continuum flow 
regime. Although the thermo chemical models em- 
ployed in continuum and particle methods are quite 
different, under conditions of thermochemical equilib- 
rium they are expected to provide identical solutions. 
The relationship between the continuum and particle 
simulation under conditions of thermochemical non- 
equilibrium, however, has not been investigated previ- 
ously. It is therefore the purpose of the present paper 
to study this relationship by computing typical hy- 
personic flows with both the continuum and particle 
simulation methods. 

Evaluation of the thermochemical models is made 
through the computation of four different cases. The 
flow conditions in the studies are given in Table 1 
and are chosen to examine the effects of vibrational 
relaxation, dissociation, and ionization. These pro- 
cesses are considered in an accumulative sense through 
a gradual increase in the initial enthalpy of the flow. 
The continuum and particle approaches employed in 
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this work are briefly described below. 


Table 1. Flow conditions. 


Case 

Ui(m/s) 

Pi (kg/m 3 ) 

Pi (tow) 

U 3 (m/s) 

1 

4000 

1.75xl0 _s 

1.17 

541 

2 

4800 

4.67xl0 _s 

31.2 

480 

3 

7310 

7.48xl0 -3 

5.00 

496 

4 

10000 

5.00xl0 -4 

0.33 

640 


Continuum Approach 

In the continuum formulation, the nonequihbrium 
gas model for air consists of eleven chemical species, 
(tf 2f 0 3> NO , N t O, N+ y 0+, iV r+ , 0 + , e~), 

and the thermal state of the gas is described by three 
temperatures: translational, rotational and vibra- 

tional (vibrational-electronic). The governing Navier- 
Stokes equations are supplemented by the equations 
accounting for thermo chemical nonequihbrium pro- 
cesses. The equation set consists of fifteen partial 
differential equations: eleven mass conservation equa- 
tions for species, one momentum equation for quasi 
one-dimensional flow, and three energy equations. For 
this work, only nitrogen reactions are considered. 
The thermochemistry model is basically that pro- 
posed by Park. 1-3 The relaxation time for vibrational- 
translational energy exchange is taken from Millikan 
and White 8 with Park’s modification which accounts 
for the limiting cross-section at high temperatures. 
Another of Park’s modifications concerning the diffu- 
sive nature of vibrational relaxation is not included 
to be consistent with the current particle model. For 
Vibration-Dissociation coupling, the average removal 
of vibrational energy due to dissociation is taken as 30 
percent of the dissociation energy. 1 The chemical re- 
action rates are prescribed by Park’s model where the 
basic dissociation rate is assumed to be governed by 
the geometric average of translational and vibrational 
temperatures. 

The numerical approach to solve the governing 
equations is fully implicit for fluid dynamics and chem- 
istry. It uses flux vector splitting for convective fluxes, 
and shock capturing. An adaptive grid strategy is also 
implemented. For the computations in this paper, a 
quasi one-dimensional code is used and a frcestrcam of 
pure nitrogen is prescribed. The details of the numer- 
ical method can be found in Refs. 4-6. 


Particle Approach 

The particle simulation code employed in this in- 
vestigation provides modeling of the translational, ro- 
tational, vibrational, and electron kinetic energy dis- 
tributions. These are complemented through simula- 
tion of dissociative, recombinative, ionising, and ex- 
change reactions. The code is vectorised for efficient 
execution on a Cray-YMP. Description of the vector- 
ised implementation may be found in Refs. 7 and 
8. The boundary conditions employed in the one- 
dimensional flow axe reflecting pistons set to the up- 
stream and downstream velocities. The downstream 
velocity is obtained either from the continuum cal- 
culations, or from available experimented data. Pre- 
vious DSMC simulations of the relaxation sone be- 
hind strong shocks 8 were commenced at the point of 
translational-rotational equilibrium (where the ratio of 
local to fceestream density is equal to 6). The present 
simulations represent an improvement in that the com- 
putation of the entire shock structure is included from 
upstream to downstream conditions. Once the shock 
reaches a specified location, small adjustments are 
made to the coordinate system of the computational 
grid to maintain a steady shock position. 

The rate of energy exchange between the transla- 
tional and rotational energy modes is simulated us- 
ing the model of Boyd. 0 The mechanics of rota- 
tional energy exchange is performed by the Borgnakke- 
Larsen 10 approach. The rate of energy exchange in- 
volving the vibrational energy mode is simulated using 
a high-temperature model also proposed by Boyd. 11 
The mechanics of vibrational energy exchange are 
computed using two different schemes. The first again 
uses the Borgnakke-Larsen approach in with a con- 
tinuous vibrational energy distribution described by a 
number of vibrational degrees of freedom, which is 
fixed. The second approach, due to McDonald, 12 al- 
lows sampling of post-collision vibrational energy lev- 
els from the discrete form of the Simple Harmonic Os- 
cillator (SHO). This approach does not require the 
value of to be estimated for the whole flowfield. 
Instead, it effectively varies ( v according to the local 
energy content of the flow, and is the preferred ap- 
proach from a physical standpoint. The manner in 
which the mechanics of energy exchange is performed 
in the particle simulation is shown by Lumpkin ef al . 18 
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Table 2. Leading constants in chemical rate data (m s /molecule /s). 


Number 

Reaction 

Continuum 3 

Particle 10 

Particle (present) 

1. 

N 3 + N 3 — > N + N + N, 

1.10x10-® T- 1 ' 8 

0.17x10"° T -1 ' 8 

7.97xl0- 18 T-°' 6 

2. 

Nj + N— *N + N + N 

4.98x10"® T -1 * 8 

t 1.85x10-* T -1 ' 8 

7.14x10-® T" 1 - 8 

3a. 

N+ + N, - N + Nj + 

1.00X1O- 1 * T°-* 

1.07xlO -17 T" 018 

1.00xlO- 18 T o s 

3b. 

N + N 3 + — ► N + + N 3 

See Ref. 1 

2.37xl0 -18 T-°-* 3 

2.34xl0- 14 T" 0 - 81 

4a. 

N + N — ► Nj + + e“ 

7.31x10-” T 11 

2.98xlO- 30 T 077 

7.31x10"” T 1S 

4b. 

Nj + + r -t N + N 

See Ref. 1 

8.88xl0 _1 ° T" 1 -” 

1.57xl0 -11 T°- 85 

5. 

N + e~ -+ N+ + e“ + e~ 

4.15x10® T-*-® 3 

l.OOxlO -14 

5.81x10-* T“ l -° 


to affect the rate of relaxation. Therefore, all the rota- 
tional and vibrational relaxation models employed in 
the particle simulations are adjusted to match the con- 
tinuum values by the correction developed in Ref. 13. 

Dissociation reactions axe modeled with the Vi- 
brationally Favored Dissociation model (YFD) pro- 
posed by Haas and Boyd. 14 As its name suggests, this 
model includes the important physical phenomenon of 
vibration-dissociation coupling. The model contains a 
free parameter <f> which controls the degree of coupling 
between vibrational and dissociative relaxation effects. 
It was demonstrated in Ref. 14 that, by increasing the 
value of it is possible to increase the dissociation 
incubation time in the simulation. Also in Ref. 14, 
through comparison with experimental data, the value 
of <j> for nitrogen was determined assuming Borgnakke- 
Larsen mechanics for vibrational energy exchange with 
a fixed value for In the DSMC code, the model 
employed for the reverse recombination reaction ap- 
propriate to VFD is that developed by Boyd. 8 All 
other chemical reactions (ionisation and exchange re- 
actions) are simulated using the steric factor developed 
by Bird. 15 The inclusion of electrons in the simulation 
is discussed in detail in Ref. 16. 

Chemical Rate Coefficients 

The rate coefficients employed in the reactions of in- 
terest in the present study are given in Table 2. These 
Me described in the usual Arrhenius form: 

k{T) = aT l exp(-E a /kT) 

where a and b Me empirically determined constants, 
E a is the activation energy, and T is the controlling 
temperature. Three different sets of coefficients are 
given corresponding to those used in the continuum 
code, in previous DSMC investigations, and in the 


present DSMC code. The values of the activation en- 
ergy used in the three sets of rate data are unchanged 
for each sepMate reaction. Therefore, the exponen- 
tial term in the Arrhenius form has been omitted from 
Table 2. 

The rate expressions employed in the continuum 
code Me those recommended in the review by Park 
et aL 3 Generally, only the forward rate constants are 
specified. In the dissociation reactions, numbers 1 and 
2, the controlling temperature in the continuum two- 
temperature approach is given by T a :=(TT„)£. For 
nitrogen dissociation, the pMticle code employs the 
rates of Byron 17 in the Vibrationally Favored Dissoci- 
ation model (YFD). It was shown previously by Boyd 8 
that these rates, when used with the VFD model, are 
capable of reproducing vibration-dissociation coupling 
observed at high temperatures. 

The reverse rates for each reaction axe obtained by 
evaluating the following temperature-dependent form 
for the equilibrium constants proposed by Park 1 : 

ln{ < K t {T y j) = A\z + Aj -f As/n(z) + A+j z + A*/ z 3 

where the A; are constants and z=10000/T. In the 
continuum code, the values of A{ for reaction 3 are ob- 
tained from Ref. 18 while those for reaction 4 are taken 
from Ref. 1. Unfortunately, this form for the equi- 
librium constant is not mathematically convenient for 
implementation in the DSMC chemistry models. How- 
ever, a set of reverse reaction rates for use in DSMC 
has been proposed by Bird, 19 and these have been used 
in a number of studies. It is possible to determine the 
equilibrium constants employed in Bird’s reaction set 
by considering the ratio of the forwMd and reverse 
ratesfor each reaction. This is performed for reactions 
3 and 4 of Table 2. The equilibrium constants em- 
ployed by Bird and those used in the continuum code 
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are shown as a function of temperature for these re- 
actions in Figs. 1 and 2, respectively. It should be 
noted that the exponential term has again been omit- 
ted for the sake of simplicity. For the charge exchange 
reaction, it is found that the equilibrium constant em- 
ployed by Bird is about 2 orders of magnitude higher 
than the continuum expression. In the associative ion- 
isation reaction, the equilibrium constant of Bird gives 
values which are again generally higher than the con- 
tinuum model. 

The goal of the present study is to evaluate differ- 
ences in the chemical models employed in the contin- 
uum and particle methods. To limit the number of fac- 
tors involved in our comparisons, it is the aim to main- 
tain consistency between the relaxation rates employed 
in the solution techniques. Therefore, a form for the 
equilibrium constant which takes the traditional Ar- 
rhenius form is fit as a function of temperature to 
Park’s expression. The limitation on the type of Ar- 
rhenius form which may be used conveniently in Bird’s 
expression for the probability of chemical reaction 1B is 
discussed by Boyd and Stark. 30 The curve fits for re- 
actions 3 and 4 are also shown in Figs. 1 and 2. The 
resulting rate constants for the reverse reactions are 
listed in Table 2. Generally, good agreement is ob- 
tained between the new DSMC expressions and Park’s 
results, particularly over the temperature range of in- 
terest, i.e. from 10,000 to 20,000 K. 

For reaction 5, the temperature dependent form 



Temperature (K) xIO ' 4 


Fig. 1 . Variation of the equilibrium constant with 
temperature for reaction 3. 


proposed in Ref. 2 is not convenient for use in the 
DSMC chemistry models. In Fig. 3, the forward rate 
constants employed by Park and Bird for this reaction 
are shown as a function of temperature. It is noted 
that Bird’s reaction rates are several orders of magni- 
tude lower. Once again, a fit is made to Park’s expres- 
sion in an Arrhenius form which may be employed in 
the particle chemistry models. The new form, which 
is given in Table 2 and included in Fig. 3, gives closer 
correspondence to Park’s results over the temperature 
range of interest. 



0.0 0.5 1.0 1.5 2.0 

Temperature (K) xIO ' 4 

Fig. 2. Variation of the equilibrium constant with 
temperature for reaction 4. 



0.0 0.5 1.0 1.5 2.0 

Temperature (K) xIO ' 4 


Fig. 3, Variation of the forward rate constant with 
temperature for reaction 5. 
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Further analyses have been performed which improve 
the correspondence between the chemical rates em- 
ployed in continuum and particle simulation for all 
reactions in air involving charged species and are re- 
ported in Ref. 16. 

Presentation of Results 


Computations are performed for four different sets 
of flow conditions for pure nitrogen and these are listed 
in Table 1 in which subscripts 1 and 2 indicate up- 
stream and downstream conditions, respectively. The 
upstream temperature is prescribed to be 300 K for all 
cases. The upstream density together with the length 
of the computational domain simulated are chosen 
such that the flows are in the near- continuum regime. 
The different upstream flow conditions also provide 
increasing enthalpy: thus, the flow behind the shock 
is characterised in Case 1 by vibrational relaxation 
processes; in Case 2 by weak dissociation; in Case 3 
by strong dissociation; and in Case 4 by weak ion- 
isation. The conditions in Cases 2 and 3 correspond 
to those investigated experimentally by Kewley and 
Hornung. 31 The results for each of these investigations 
are described in the following sub-sections. 


Case 1: Vlbrationally Relaxing Flow 

Density profiles for the first case investigated are 
shown in Fig. 4. Very good agreement is found be- 
tween the continuum and particle siumulation results. 
Two different DSMC computations are shown: the 
first corresponds to the use of the B or gnakke- Larsen 
approach (BL) for performing the mechanics of vibra- 
tional energy exchange with a constant number of vi- 
brational degrees of freedom, £*=1.0. This value cor- 
responds closely to that evaluated at the downstream 
equilibrium temperature. The second solution em- 
ployed the discrete vibrational energy sampling ap- 
proach for the Simple Harmonic Oscillator (SHO) of 
McDonald 7 which automatically varies This is the 
first time that comparison is made between continuum 
and particle simulations for vibrational relaxation be- 
hind a strong shock. It is very encouraging to observe, 
under near-continuum conditions, that the two meth- 
ods give such close agreement. 


fS 
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Fig. 4. Comparison of continuum and particle solu- 
tions of the local to upstream density ratio 
for Case 1, 
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Fig. 5. Comparison of continuum and particle so- 
lutions of the translational and vibrational 
temperatures for Case 1. 


The variation in translational and vibrational tem- 
peratures for this case are shown in Fig. 5. The parti- 
cle solutions are obtained with McDonald’s variable ( v 
model. Once again, very good agreement is obtained 
between the two solution techniques. Temperature is 
generally a much more sensitive quantity to simulate 
than density. The close correspondence between the 
continuum and particle results indicates that the vi- 
brational relaxation models employed in each are very 
nearly equivalent. This comparison therefore lends 
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strong support to the use in the particle simulation 
of the vibrational energy exchange probabilities devel- 
oped for the DSMC method, 11 the correction term 
required to equate the continuum and particle relax- 
ation rates, 13 and the mechanics of vibrational energy 
exchange. 13 It should be noted that the degree of dis- 
sociation under these flow conditions is less than 1%. 


Case 2i Weakly Dissociating Flow 

The second set of conditions considered has an in- 
creased enthalpy which gives rise to weak dissociation 
behind the shock. This case is of additional inter- 
est as it was studied experimentally by Kewley and 
Hornung 31 who employed interfero grams to measure 
the variation in density behind strong shocks of nitro- 
gen. The increase in enthalpy is revealed in the density 
profiles shown in Fig. 6 in which the normalised den- 
sity rise reaches a value of about 10 at the downstream 
boundary. While both solutions give good agreement 
with the experimental data, it is clear that the parti- 
cle solution provides the better correspondence. The 
DSMC profile is obtained with the variable model. 
Comparison of the translational and vibrational tem- 
peratures computed through the shock are shown in 
Fig. 7. Again, a very good agreement is observed for 
the two sets of results. Considering the excellent agree- 
ment obtained in Fig. 7 between the continuum and 
particle methods, and also for the case of vibrational 
relaxation, it is concluded that the differences observed 
in Fig. 8 must be due to the dissociation models em- 
ployed in each simulation technique. This indicates 
that the continuum two- temperature model gives a 
dissociation rate which is slightly slower than that of 
experiment and DSMC, In other words, for a weakly 
dissociating gas, while the effect of dissociation on vi- 
brational relaxation is small, the effect of vibrational 
relaxation on dissociation is overestimated in Park’s 
two-temperature model. 

The results for the mole fractions of molecular and 
atomic nitrogen are shown in Fig. 8. As expected from 
the previous comparisons, there is close agreement be- 
tween the two numerical approaches with DSMC pre- 
dicting slightly more dissociation than is obtained in 
the continuum solution. 
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Fig. 6. Comparison of continuum and particle solu- 
tions of the local to upstream density ratio 
for Case 2. 
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Fig. 7. Comparison of continuum and particle so- 
lutions of the translational and vibrational 


temperatures for Case 2. 


Case 3: Strongly Dissociating Flow 

The further increase in enthalpy for Case 3 gives 
rise to stronger dissociation effects. Once again, the 
flow conditions modeled match those considered exper- 
imentally by Kewley and Hornung. 31 The experimen- 
tal profile of density behind the shock is compared with 
several different computational results in Fig. 9. The 
comparison between the continuum solution and the 
experimental data is excellent. Two different DSMC 
profiles are also shown in Fig. 9: the first employed 
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Fig, 8. Comparison of continuum and particle solu- 
tions of the species mole fractions for Case 2. 
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Fig. 9. Local to upstream density ratio for Case 3: 
particle simulation employed old YFD and 
exchange mechanics models. 


and ^=3. This represents the model proposed 
previously in Ref. 8. The second particle simulation 
employed the variable model and again used <p=Z. 
For both of these DSMC modeling configurations, the 
computed density profile lies a little below both the 
experimental data and the continuum solution. A fur- 
ther DSMC solution is generated in which the YFD 
parameter (p is reduced to 2 for use with variable £ c . 
The results of this simulation are compared with both 
the experimental and continuum profiles in Fig. 10. As 
expected, a reduction in the value of <p reduces the 



Fig. 10. Local to upstream density ratio for Case 3: 
particle simulation employed new YFD and 
exchange mechanics models. 


degree of vibration- dissociation coupling which leads 
to a slightly increased dissociation rate immediately 
behind the shock. With this model configuration, 
the particle method provides excellent agreement with 
both the experiment and the continuum solution. It 
should be noted that this is the DSMC model config- 
uration employed in the computations for Case 2. 

The translational and vibrational temperature pro- 
files computed with the continuum and DSMC tech- 
niques are shown in Fig. 11. Generally, very good 
agreement is observed between the two. There is a 
noticeable difference in the peak vibrational tempera- 
tures computed by the two methods. This has quite 
significant implications for the estimation of radiative 
emission in such flows. The difference is attributable 
to dissociation- vibration coupling, i.e. how the vibra- 
tional energy distribution is affected by dissociation. 
This process is modeled quite differently in the con- 
tinuum and particle approaches. These results indi- 
cate the need for experimental measurement of vibra- 
tional temperature profiles behind shock waves under 
conditions similar to those considered here. For com- 
pleteness, the profiles of mole fractions of the neutral 
species are shown in Fig. 12. The stronger degree of 
dissociation for these flow conditions is very evident 
and, as expected, very good agreement is found be- 
tween the two solutions. 
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Pig. 12. Comparison of continuum and particle solu- 
tions of the species mole fractions for Case 3. 


For this strongly dissociating case, it is found that 
Park’s two temperature model reproduces the experi- 
mental data very well. It is very encouraging that the 
two temperature model gives such a favorable com- 
parison with the experimental data in strongly disso- 
ciating flow as this is the regime for which the model 
has been developed. Indeed, the present comparison 
arguably provides the strongest evidence to date that, 
despite its weak theoretical basis, the two-temperature 
model does produce adequate simulation of strongly 
coupled vibration-dissociation processes. The present 


investigation is unique in that evaluation of the model 
is performed through direct comparison with experi- 
mental measurements of a fundamental flow quantity. 
The mqdel was previously calibrated against experi- 
mental data by Park 1 through comparison with ra- 
diation emission spectra, and by Candler 33 through 
comparison with shock stand-off distance. Due to the 
excellent comparisons between DSMC and experiment 
in Figs. 0 and 10, it is recommended that McDonald’s 
collision mechanics and the VFD model with <f )~ 2 be 
used for simulating nitrogen dissociation with the par- 
ticle method. 

Case 4; Weakly Ionising Flow 

The increase in enthalpy for Case 4 is sufficient to 
give rise to significant ionization effects behind the 
shock. In performing the DSMC computations of 
the ionized flowfield, a steady shock solution is first 
obtained with the ionizing reactions omitted. After 
reaching this point, the ionized species are included 
and a further short transient phase in the simulation 
then allowed before sampling of flow properties is com- 
menced. These procedures are adopted because the in- 
clusion of electrons in the flowfield requires a reduction 
in computational time-step by two orders of magni- 
tude. To compute the entire flowfield with such a small 
time-step would expend much larger computational re- 
sources. The comparison for density profiles computed 
with the numerical techniques is shown in Fig. 13. As 
with the previous cases, good agreement is obtained 
between the two solutions. The temperature profiles 
computed with the continuum and particle methods 
for the translational and vibrational modes axe com- 
pared in Fig. 14. The peak values for each energy mode 
are in good agreement. It is observed that the trans- 
lational temperature shock computed with DSMC is 
thicker than the continuum result. This is due to the 
relatively low upstream density employed in this inves- 
tigation. A more thorough analysis of such behavior 
will form the basis of future study. The computed mole 
fractions for the neutral species obtained with the nu- 
merical techniques are compared in Fig. 15 and those 
for the charges species are compared in Fig. 10. The 
agreement which is generally obtained is very satisfac- 
tory. This comparison verifies that the new forms of 
the reverse reaction rates employed in the particle sim- 
ulation are nearly equivalent to the expressions used 
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Fig* 14. Comparison of continuum and particle so- 
lutions of the translational and vibrational 
temperatures for Case 4. 
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Fig. 15. Comparison of continuum and particle solu- 
tions of the neutral species mole fractions for 
Case 4. 



0.0 0.5 1.0 1.5 2.0 2.5 


Distance (cm) 

Fig. 16. Charged species mole fractions for Case 4: 
particle simulation employed new chemical 
rate data. 


in the continuum analysis. It should be noted that a 
degree of statistical scatter is exhibited by the DSMC 
results for the less abundant species. To assess the 
effect of using the new reaction rates > a particle simu- 
lation is also performed with Bird’s rate data 15 . The 
variation in the mole fractions of the charged species 
computed in this way are compared with the contin- 
uum results in Fig. 17. None of the species profiles arc 
found to be in good agreement. With Bird’s rates, the 
most populous ion is N* + , whereas the new particle 


rate data agrees with the continuum solution in giv- 
ing N + as the most abundant ion. With Bird’s rate 
data, the mole fraction of electrons at the downstream 
boundary is about 2.5xl0 -3 whereas the continuum 
simulation gives a value of about 1.8x10 2 . If it is 
accepted that the rate coefficients provided in Refs. 1 
and 2 are the more physically realistic, these large dif- 
ferences observed with Bird’s older data set must call 
into question previous DSMC investigations which em- 
ployed those reaction rates. 
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Fig. 17. Charged species mole fractions for Case 4: 
particle simulation employed old chemical 
rate data. 

Concluding Remarks 

This study was motivated by the requirement to 
evaluate the relationship between continuum and par- 
ticle simulations of hypersonic flows in the near- 
continuum regime. The results obtained in the in- 
vestigation have established that a close correspon- 
dence exists between the thermochemical nonequilib- 
rium models employed in these solution techniques. In 
the case of vibrational nonequilibrium, the agreement 
between the two sets of numerical results validated a 
number of recent modeling developments for comput- 
ing the rate and mechanics of vibrational energy ex- 
change in the particle simulation. In the cases of weak 
and strong dissociation, both the continuum and par- 
ticle models for vibration-dissociation coupling were 
successfully evaluated against experimental data. This 
is a most interesting result considering the large differ- 
ences in the dissociation models employed in the two 
techniques. In the case of weakly ionised flow, it was 
necessary to develop new forms for some of the chem- 
ical rate constants for use in the particle simulation. 
These were developed so as to be nearly consistent 
with the continuum expressions, and also to be math- 
ematically convenient for use in the particle chemistry 
models. The next stage in this continuing investigation 
will be evaluation of these methods for flow conditions 
in the transition regime, i.e. at higher Knudsen num- 
bers. In such flows, rarefaction effects may invalidate 


use of the Navier-Stokes equations, and give rise to 
large differences between the continuum and particle 
simulation results. 
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Abstract: 

This paper describes the computational work be- 
ing performed at Ames on the simulation of the 16” 
Shock Tunnel facility. The paper describes the ap- 
proach used and shows some preliminary results for 
various flow transients. In particular, we describe the 
numerical problems encountered during the compu- 
tation of these flows, and the methods used to resolve 
them. We also discuss the validity of some approx- 
imations used, notably concerning the reduction of 
the problem into problems of smaller dimensionality, 
or smaller size. We show how quasi-ID simulations 
can be used to help design experiments, or to better 
understanding the characteristics of the facility. An 
application to the design of a non-intrusive diagnos- 
tic is shown. The multi-dimensional flow transients 
computed include the shock reflection at the end of 
the driven tube, the shock propagation down the noz- 
zle, and the breaking of the main diaphragm. The 
interaction between separate flow events will also be 
discussed. 


I. Introduction 

The Ames 16” facility, shown schematically in 
Figure 1, can be considered as a typical example of a 
shock- tunnel for hypersonic flows. The shock tunnel 
is about 70 meters long, composed of a driver sec- 
tion (17” diameter) filled with a light combustible 
mixture, and a long driven tube (filled with the test 
gas at low pressure). At the end of the tunnel is a 
supersonic nozzle 6 meters long (1 m diameter at the 
exit plane) and finally a test section. The complete 
operation of the facility typically results in test times 
of the order of 5 to 20 milliseconds: the time for the 
main shock to propagate from the main diaphragm 
to the end of the driven tube is of the order of 7 msec. 
After partial reflection at the end of the driven tube, 
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the shock interacts with the incoming contact discon- 
tinuity (CD) which separates the test gas from the 
driver gas. For ’tailored’ conditions, there is no wave 
reflected back from this interaction, and the flow be- 
tween the CD and the nozzle is uniform and steady. 
More detailed descriptions of a shock tunnel opera- 
tion can be found for example in ref. [1]. The test 
time is measured from the time the flow conditions 
in the stagnation region become steady, until the CD 
reaches the end of the driven tube, and the nozzle 
flow becomes contaminated with the driver gas. A 
flow ’quality’ will be represented by the steadiness of 
the stagnated flow, the low level of contamination of 
the test gas by driver gas or impurities, as well as 
the peak conditions (pressure, enthalpy) attainable 
by the facility. 

Numerical simulations are required to better un- 
derstand the test flow conditions (temperature, mass 
fraction profiles, time dependence, etc..) and sup- 
plement the experimental knowledge: they are used 
also to improve the design of the experiments or to 
improve the tunnel operation characteristics. In or- 
der to satisfy these objectives, several aspects of the 
tunnel operation need to be simulated. Modeling 
the entire facility from the driver to the test section 
requires us to grid a physical length of about 50 me- 
ters; yet, important flow features such as shock and 
contact discontinuities (CD) should be resolved with 
good accuracy, i.e. a few mm in the axial direction. 
The spatial stiffness is then of the order of 10 4 ; this 
is a conservative estimate, since the boundary layers 
also need to be resolved in some cases, requiring grid 
spacing as low as tens of micrometers in the radial 
direction. The modeling of the shock-tunnel facil- 
ity falls into the general category of multiple scale 
problems, and calls for appropriate special numeri- 
cal methods which we will hint at in this paper. 

The complete operation of the hypersonic facility 
also involves a large number of different physical pro- 
cesses at work, some of which are not well understood 
or are very difficult to compute. The combustion 
process in the driver gas requires a good description 
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of the energy deposition and flame propagation; the 
main diaphragm rupture would require us to under- 
stand the material deformation up to the plasticity 
limit. The penetration of the'jet of driver gas into 
the driven tube is a problem of 3D, turbulent, multi- 
scale mixing. The temperatures are sufficiently high 
that some wall ablation takes place and contaminates 
the flow. Finally, radiative effects may need to be 
considered, and low density, thermal nonequilibrium 
gas models may be required in the nozzle expansion. 
A multiplicity of CFD tools must therefore be made 
available to study the effects of each of these phe- 
nomena. 

The problem is compounded by the lack of cru- 
cial experimental data: such phenomena such as di- 
aphragm rupture are difficult to observe, and only 
rough estimates of the process time and energy scales 
can be made available and approximately correlated 
with the experimental data. Ablation of the tun- 
nel material is also a very complex physical process, 
which may depend on the microstructure of the ma- 
terial itself. 

Solving this problem will require us to develop 
and test appropriate numerical techniques, decom- 
pose and reduce the problem accordingly, and de- 
velop and test various physical models. Simplifica- 
tion of the problem can be accomplished first by ’di- 
mensional reduction’, i.e. solving the problem in in 
a reduced number of dimensions. Quasi-ID simula- 
tions have low computational time and memory re- 
quirements, and therefore can be used for a large 
number of computations, such as those required for 
design and sensitivity studies. It will be important 
to verify the accuracy of these simulations, and to 
estimate the quality of information which can be ob- 
tained from them. 

In a typical cartesian tradition, decomposition 
of the problem into several independent sections is 
usually attempted. This we regard as ’causality re- 
duction’. This approach requires several approxima- 
tions and simplifications, and its validity needs to be 
more clearly justified. For example, the time evolu- 
tion of the flow in the test section clearly requires 
us to accurately compute the unsteady flow in the 
stagnation region, including gas/wall interactions, 
shock/boundary layer interactions. The shock/CD 
interaction is also very important, and this requires 
us to know very well the shape of the CD. Going 
back further, the shape of the CD will be influenced 
by the earlier process of flow establishment from the 
main diaphragm rupture and thereon. But the de- 


tails of the diaphragm rupture will depend on the 
pressure rise in the driver gas, and this requires us 
to model first the turbulent flame propagation in the 
driver tube. The logical sequence of events, i.e. from 
driver tube combustion to flow expansion into the 
test chamber, is unfortunately the most difficult to 
compute, since it requires us to be able to model the 
most complex phenomena at first. Therefore we have 
to sacrifice some degree of detail in order to obtain 
practical answers in reasonable time. By causally re- 
ducing the problem, we loose accuracy but gain in 
efficiency: the method can still be used to examine 
the important physical effects (i.e. sensitivity stud- 
ies) in an efficient way. 

In this paper, we will discuss these methods and 
the preliminary results. The following questions 
must be addressed: 

• what improvements in numerical techniques are 
required. 

• how accurate are the results from reduced (ID) 
models. 

• how significant are the causal influences of var- 
ious sections of the tunnel. 

« what physical phenomena are important and 
how to model them. 

The last item on our list will not be discussed 
here. In this paper, we will concern ourselves primar- 
ily with the numerical techniques, and the evalua- 
tion of the approximations commonly made in shock- 
tunnel simulation. The obvious intent of the devel- 
opment of this numerical capability is to provide an 
array of numerical tools for better understanding of 
the facility operation, and the design of new test con- 
ditions and diagnostic procedures. The combined fo- 
cus of both experimental and numerical tools leads to 
superior measurement capability, and is an approach 
practiced at NASA- Ames [2]. 

II. Numerical Techniques 

We have already mentioned that the spatial and 
time stiffness can be quite large; fortunately the CFD 
community has developed an arsenal of techniques 
which can potentially be used to reduce the severity 
of the problem. First, we observe that the flow dis- 
continuities are few and generally well localized: it 
is therefore unnecessary to simultaneously carry the 
same resolution requirements in all regions of the fa- 
cility. We can rely on several techniques of dynamical 
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grid refinement to carry this task. Second, the prob- 
lem can be approximately separated in three parts. 
The combustion process in the driver gas proceeds 
independently of the driven tube, and can be com- 
puted separately if necessary (throughout this paper, 
it will be simply modeled). Since the flow is mostly 
supersonic in the nozzle, there is no backward influ- 
ence of the flow from the nozzle to the tubes. It is 
therefore adequate to compute at first the unsteady 
flow in the driver and driven sections, including the 
nozzle and up to the throat; the flow solution at the 
exit plane can be stored according to the required 
spatial and time accuracy, and used as unsteady in- 
flow condition for the remainder of the facility, i.e. 
nozzle and test section. This allows us also to switch 
between different physical models (e.g. nonequilib- 
rium thermodynamics) and different dimensionalities 
(3D versus 2D axi-symmetric or 2D versus ID). 

Another important question concerns the choice 
of numerical method: although we are interested 

in the unsteady flow motion, the flow time scales 
that need to be resolved may be large compared to 
the time scale from the stability limit of an explicit 
scheme. This is true especially for a clustered vis- 
cous grid. An implicit method, if time-accurate, can 
potentially be used for this problem as well. How- 
ever, there is a huge cost associated with the inver- 
sion of the block-tridiagonal matrices from the LHS 
of the system of equations, especially when multi- 
ple species are present. Our experience showed that 
even when only two species are considered (driver 
gas/driven gas), the implicit method needs to be ex- 
ecuted at a CFL number > 5 for better performance 
than the explicit method. Due to the severe nature of 
the problem (pressure ratio, shock speed), this CFL 
value is generally too large for these transient flows. 
However, the viscous terms alone can be treated im- 
plicitely, and show greater stability than the convec- 
tive terms: our method therefore performs an Opera- 
tor Splitting between convection, viscous dissipation 
and chemistry. The viscous terms are treated im- 
plicitely, and the inversion of the 3 x 3 (in 2D) block 
tridiagonal matrices is sufficiently fast to be justified. 
The chemistry is also treated point-implicitely, and 
the algorithm is also very fast. When the chemistry 
is extremely stiff however, linearization errors will 
reduce the stability of the point-implicit algorithm. 
For example, let us consider air (79% No, 21% Oo) 
suddenly raised to 100 atm and 6000 °K: at con- 
stant volume, the system reaches equilibrium within 
a fraction of a microsecond. Attempting to solve 
implicitely the chemistry in one iteration and with a 
time step larger than ~ 2 10~ 5 sec would lead to very 


unstable and unphysical solutions. The initial evo- 
lution from the highly nonequilibrium state needs to 
be time-resolved more accurately: our algorithm sub- 
iterates (more precisely, ’sub-cycles’) the chemistry, 
i.e. integrates over a given time interval At with vari- 
able steps St. The sub-step is initially estimated from 
the rate of change of the species, then is stretched in 
a fixed proportion for the subsequent sub-iterations. 

1 - (1 + a)St n , with a > 0. Each sub-iteration 
is solved point-implicitely. The change in tempera- 
ture is also estimated at each sub-iteration from the 
change in chemical energy. Figure 2 shows a sample 
computation (a = 0.5) for this highly nonequilib- 
rium case. The ’exact’ time evolution is also shown. 
The integration intervals are At = 5 10~ 7 seconds. 
One can see that the variable step solution is off in 
the equilibrium region until the end of the first inte- 
gration interval, at which time the thermodynamic 
state of the gas is re-analyzed and a more accurate 
calculation of the temperature is performed. The 
sub-iterated solution then quickly matches the exact 
solution. 

To maintain good accuracy, the relative changes 
in temperature estimated during the chemistry 
should be limited. A chemical time scale is there- 
fore defined as the time during which the chemical 
reactions induce a relative change in temperature of 
~ 5%f . The effect of chemistry on the flow dy- 
namics occurs principally through the temperature 
change, and this time step limit provides us with 
a criterion for accurate coupling of the chemistry 
and flow dynamics. If there is no monitoring of the 
chemical time scale and no enforcement of this time 
step limit, an instability generally develops rapidly, 
especially in cases of stiff chemistry; this instabil- 
ity can occur for both a fully-coupled approach or 
the Operator-Splitting approach. The case of very 
stiff chemistry would seem therefore to be very in- 
efficient: this corresponds basically to a very poor 
resolution of the chemical relaxation distance, such 
as the one behind the primary shock. However, it 
is possible to somewhat rescale the chemical time 
scale without changing the solution too much. Let 
us suppose that we want to resolve the flow dynamics 
on the order of a convective time-scale At conv , but 
the chemistry is so stiff that accuracy and stability 
considerations restrict us to (for example) a global 
time scale — 0.01A2 CO nv We can artificially rescale 
the chemistry by 10, and the stability limits lead us 

^This number may be varied, according to the desired 
accuracy and/or the stability. It is our experience that in 
most cases, the maximum relative change allowed should be 
less than 10%. 
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to a time scale 0.1At conv , which we can better af- 
ford. This rescaling can be done locally, depending 
on the stiffness ratio: this will affect slightly the so- 
lution, for example by over-estimating the chemical 
relaxation length behind a shock, but the difference 
may be practically insignificant for very stiff cases. 
The technique is easily implemented by restricting 
the integration interval (it becomes 5 10“ 8 sec in our 
example): the insert in Figure 1 shows the results 
of this rescaling approach. We see that the original 
transients are still reproduced, and the solution still 
converge to the exact solution. This sub-iterated al- 
gorithm, with or without the rescaling option, is used 
in all our calculations. 

The spatial stiffness can be solved by various tech- 
niques of dynamic grid adaption. In the driver tube, 
strong gradients and discontinuities exist only dur- 
ing the combustion process. After rupture of the 
main diaphragm, only weak gradients subside, and 
this section of the facility does not require high reso- 
lution. The spatial accuracy is required for the prop- 
agating shock and CD, down the driven tube. One 
possible technique is to dynamically compress and 
stretch the grid points to accumulate them in the 
required regions. This technique can work well in 
one dimension and is easy to implement: care must 
be taken however to prevent sudden jumps in spac- 
ings, or accuracy will be lost. The technique may be- 
come more problematic in two or more dimensions, 
as some flow features (including the CD) may have 
more complex shapes: the grid can become distorted, 
and can affect the solution. Another technique, pio- 
neered by a group at Livermore [3], consists of adding 
smaller scale grids in regions of interest. These grids 
can be exact sub-scale replicas of some regions of 
the coarse, background grid: the flow variables can 
be transferred between grids in an exactly conser- 
vative way. This way, the grids are never distorted 
and their motion can be computed very quickly, with 
minimal overhead. An example is shown in Figure 3: 
the propagation of the primary shock and CD down 
the driven tube is computed in a one-dimensional 
model. A first calculation used a constant spacing 
grid, with about 1000 points for the whole tunnel. 
A second calculation used a background grid for the 
whole tunnel with 250 points, and a high resolution 
grid super-imposed on it: the latter moves along with 
the shock, with a sliding motion. The sub-scaled grid 
had a size of 800 grid points, and the scale ratio be- 
tween the two grids was 20 to 1. Flow values within 
the background cells are computed by volume aver- 
aging of the sub-scaled cells present, if any. While 
both calculations used approximately the same num- 


ber of points, it is clear that the second method has 
a much higher local resolution, and the sharpness of 
the CD is dramatically improved. 

Another method, which we will be testing in the 
future, uses non conservation of the number of grid 
points. Unstructured grids can be manipulated to 
accumulate points in the regions of interest, either 
by grid displacement or by dynamic subdivision. 

The singular axis was found to be another recur- 
rent problem during the computation of flow tran- 
sients: the simulation of the shock reflection at the 
end of the driven tube, for example, initially showed 
a strong conical shock structure near the axis with 
the apex of the cone leading the overall structure (see 
Figure 4). This peculiar formation can also be seen 
for example in the results of P. Jacobs [4]. This struc- 
ture is possible only if a very intense and high veloc- 
ity jet of gas is produced and maintained on the axis: 
this is a highly singular and unphysical behaviour. 
Close examination of the numerical results indeed 
showed an excessively large axial velocity component 
in a single cell close to the axis. Because this high 
velocity jet was present in one cell only, and did not 
show signs of diffusion, we were convinced that it 
is the result of a numerical error. A similar phe- 
nomenon can be observed also in the propagation of 
the main shock down the nozzle [5], and can be seen 
at the axis or in some cases near the walls. This phe- 
nomenon was observed only for axi-symmetric flows, 
and when a second-order accurate scheme (with min- 
imal dispersion) was used. It was finally related to 
the aspect ratio of the grid cells: indeed, the problem 
disappeared when the aspect ratio of the grid cells 
was adjusted to lower the radial gradients. If the 
spacing is such that Ax << Ar, the flow features 
are smooth and accurate. If Ax > Ar an instability 
may develop in some regions of the flow. We assume 
the problem comes from the axi-symmetric pressure 
correction term, which is not part of the monotonic 
(TVD) fluxes and acts as a non-conservative momen- 
tum source term. This pressure correction term is a 
result of the formulation of the Euler equations for 
an axi-symmetric problem, and can be easily visu- 
alized with the following finite-volume description. 
Let us consider a cell in a cylindrical geometry, de- 
scribed schematically in Figure 5. As one approaches 
the axis, the ratio of cell side surface to cell volume 
behaves as 1/Ax for the axial direction, while it is 
exactly 0 for the radial direction. There is a contribu- 
tion to the radial momentum density from the pres- 
sure on the sides of the cells in the azimuthal direc- 
tion. This contribution is proportional to 1/Ar. It is 
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clear that when the cells are clustered near the axis, 
the contribution from the non-conservative momen- 
tum source term may become dominant, and there- 
fore there is no guarantee that a non-oscillatory solu- 
tion can still be obtained. The problem usually dis- 
appears when the grid is relatively coarse in the ra- 
dial direction, near the axis, in which case the source 
term ceases to be the dominant contribution. This 
observation also applies to other cases, such as the 
computation of shock layers on blunt axi-symmetric 
bodies. In most unsteady cases or for steady flows, 
the radial gradients will be weak near the axis. One 
particular exception would be the propagation of a 
shock ring towards the axis: this is such the case 
in the reflection of the primary shock in the tunnel, 
from the upper section of the end wall. Figure 6 
shows a time sequence for the same case as Figure 4, 
but for a grid coarsened near the axis: we see that 
the instability is effectively removed. A similar im- 
provement can be obtained by reducing the spatial 
order of accuracy in that region. 

The general effect of the grid structure or aspect 
ratio on the flow solution should never be under- 
estimated, especially in the severe cases we are con- 
cerned with. Although we may see unstructured 
grids as a potentially useful tool for the flows of in- 
terest here, their effect on the flow solution will have 
to be carefully evaluated. 

The flow expansion in the supersonic nozzle is one 
of the easier tasks to perform: the transient flow sim- 
ulation is required to examine the steady flow estab- 
lishment time in the test section, and to verify that 
the flow does not choke during that time. This is 
particularly relevant for large objects or for massive 
gats (fuel) injection. Since the nozzle is evacuated to 
very low pressure before the rupture of the secondary 
diaphragm, the flow propagating down the nozzle is 
preceded by a high velocity jet. The mean free path 
in that region can be quite large (~ 3 mm for P = 100 
milliTorr), and strictly speaking the Navier-Stokes 
equations cease to be valid in this low-density re- 
gion prior to the shock. By attempting nevertheless 
to solve the flow dynamics with the Navier-Stokes 
equations, we are experiencing strong viscous effects 
which operate on a very short 10 pico-second) 
timescale. Although the viscous fluxes are computed 
implicitely, it is necessary for stability reasons to ar- 
tificially reduce the strength of these viscous effects 
in that region. This is easily accomplished by us- 
ing a numerical switch that effectively and smoothly 
removes the gradients on a scale smaller or compa- 
rable to the mean-free path during the calculation 


of the viscous fluxes. This switch will effectively re- 
duce the spreading of the shock into the low density 
gas. The shock thickness will therefore be under- 
estimated. The cutoff length scale can be adjusted 
arbitrarily: in Figure 7 we sho^v the effect of the cut- 
off choice (0.1 versus 10 mean free paths) on the solu- 
tion, for two cases of nozzle pressures. The change is 
dramatic for the first case (100 //Torr) and unnotice- 
able for the second (100 milliTorr). Nevertheless, for 
the latter case, which is also the case we will be using 
throughout the remainder, the viscous time scale is 
reduced by a factor of 5 by choosing the higher (10 
m.f.p.) cutoff value. 

Using all the modified techniques now at our dis- 
position, we are able to sucessfully compute the tran- 
sient flow’s in the nozzle and driven tunnel at the 
rupture of the secondary diaphragm (see Figure 6). 
Different grids were used for the driven tube and 
nozzle regions, and no subscaled grids w r ere used for 
this particular case: both regions were coupled at 
each iteration. The results of the computation for 
the nozzle flow (Figure 8) are in good qualitative 
agreement with an experimental schlieren of a simi- 
lar problem shown in Figure 9, taken from ref. [6]. 
Both show the curved leading shock, and a complex 
structure behind it, dominated by a Mach disk, it- 
self supported by oblique shocks emanating from the 
nozzle walls. Although the initial conditions for these 
two flows are very different, the similarity between 
the two structures lead us to conclude that they are 
examples of a general pattern. Details of the con- 
ditions (pressures, geometry, etc..) may change the 
relative strength and position of the shocks, with- 
out modifying the overall configuration. A snaphost 
taken at later times would show that near the exit 
plane of the nozzle, the primary shock straightens 
out and becomes normal. If one was to perform an 
unsteady computation in the test chamber assuming 
a normal shock at the inflow, the calculation would 
be in error by leaving out the complex shock struc- 
ture which propagates immediately behind the lead- 
ing shock. This is shown for example in Figure 10, 
where a cone has been used as testing body. 

III. Dimensional Reduction 

When test times are large compared to the tran- 
sients, it seems appropriate to compute the nozzle 
flow for the steady state, without having to perform 
the calculation time-accurately for tens of millisec- 
onds. This allows us to use many numerical tech- 
niques to make this computation more efficient. We 
have done several such computations, solving for the 
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full Navier-Stokes equations, but for limited regions 
at a time. The computation proceeds on a subset 
of the whole nozzle grid, until convergence is ob- 
tained. The subset then moves further down the 
nozzle; the procedure is similar to the displacement 
of a ’window’ along the flow direction. This proce- 
dure is half-way between a global calculation and the 
PNS method; although less efficient than the PNS 
method, it is more flexible and allows us to correctly 
treat embedded recirculation zones or other unex- 
pected subsonic regions which may occur. It was 
observed that the convergence rate and the compu- 
tation time are dominated by the chemistry in the 
boundary layer, and the downstream propagation of 
chemical changes in the low velocity sublayer. Al- 
though the calculation does not display severe in- 
stabilities such as may be the case for the transient 
flows, it is still desirable to coarsen the grid in the 
radial direction near the axis, to avoid spurious fluc- 
tuations in pressure. 

The steady nozzle flow computations are usually 
done to obtain estimates of the flow conditions at the 
exit plane, in particular temperature, Mach number, 
pressure and species concentrations. Since it is gen- 
erally assumed that the flow is nearly uniform at the 
exit plane, it seems that simple one-dimensional com- 
putations should be adequate. We have performed 
such computations using a quasi- ID version of our 
code, and compared the results. The effect of the 
boundary layer is expected to be the major contri- 
bution to any potential disagreements between the 
quasi- ID solution and the 2D axi- symmetric solu- 
tion. This effect is taken into account by effectively 
modifying the nozzle profile and its area. By con- 
sidering the inviscid core only, the quasi- ID method 
can approximately take into account the constricting 
effect of the boundary layer by assuming a new noz- 
zle shape for which the effective radius is a constant 
fraction of the real nozzle radius. The dependence 
of some flow quantities on the effective nozzle area is 
shown in Figure 11. It is clearly apparent that the 
static pressure has the highest sensitivity to the ef- 
fective nozzle area, and therefore to the boundary 
layer thickness. Temperature and species concen- 
trations are less sensitive to the area variation, and 
additional uncertainties about the chemical rates or 
contaminants are likely to cloud the issues: the static 
pressure is therefore an important variable to mea- 
sure at the nozzle exit. The comparison with the 2D 
axi-symmetric results is shown in Figure 12, where 
again the sensitivity of the static pressure is clearly 
demonstrated. The best agreement is obtained for 
an effective nozzle radius of 87% the real radius (i.e. 


a boundary layer thickness approximately 13% the 
nozzle radius). Notice also that the static pressure 
is the only variable that still displays significant ra- 
dial fluctuations within the inviscid core, at the exit 
plane. A multi-point measurement of this variable is 
therefore doubly informative, as far as code valida- 
tion is concerned. In Figure 13 we show the radial 
profile of Temperature and NO concentration at sta- 
tion code-named N3 (2.37 meters downstream of the 
throat) and the exit plane. All the species profiles 
are relatively unchanged between N3 and the nozzle 
exit, except atomic nitrogen N } which is present in 
very small amounts. Owing to the high reactivity of 
N, this result is not surprising: most of the chem- 
istry is frozen before station N3, and this is especially 
true of NO. The temperature variation between N3 
and the exit plane is therefore due to hydrodynami- 
cal effects only. 

NO is a relatively important specie to measure 
(it effectively ties up a significant amount of oxygen 
and has a noticeable effect on ignition delays), and it 
can be done easily with a laser absorption method. A 
computational study of this diagnostic technique can 
also be done to help design the experiment. The com- 
puted high resolution spectrum of the NO( y) band 
system is shown in Figure 14, with and without the 
boundary layer. The intensity at peak absorption of 
the (0, 0) band changes by a factor of ~ 3 when the 
boundary layer is removed. A precise measurement 
of the core flow must avoid the uncertainty caused by 
the boundary layer: this can be simply done by pro- 
viding an optical wave guide that protects the beam 
from the boundary layer. Although this protrusion 
is likely to disrupt the flow, the effect may be small 
enough or irrelevant at the exit plane. We can now 
look at the effect of core flow variation by comparing 
the intensity for the 2D solution (without BL) and 
a quasi-ID solution. The result is plotted in Fig- 
ure 15-a as the ratio of intensities between the two 
cases. Surprisingly, the error is less than 10% for the 
most sensitive lines. At the exit plane (Figure 15-b), 
the results are even better, with deviations less than 
2%! This seems to indicate that this measurement 
can yield useful information on the core flow, with- 
out having to worry too much about details of the 
radial fluctuations within the inviscid core. 

IV: Causality Reduction 

In the previous section we have seen that dimen- 
sional reduction of the problem can still yield some 
useful results, and quasi- ID calculations can help de- 
sign and understand the facility operation and exper- 
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imental diagnostic techniques. Some flow features 
are lost in this reduction process, for example the 
shock structure that propagates down the nozzle. If 
one were to compute the transient flow inside the test 
section, one would require a prior computation of the 
transient nozzle flow. This is an example of causal 
interaction between two different regions of the fa- 
cility. The computation of the transient nozzle flow 
can be done by assuming constant stagnation con- 
ditions, or starting from a uniform shock impinge- 
ment on the end wall of the driven tube. For later 
times, the propagation of the reflected shock in the 
driven tube must be followed accurately: this entails 
the interaction of the shock with the bounday layer 
and with the contact discontinuity. If this numeri- 
cal simulation is done with sufficient detail, we can 
gather information on the effect of boundary layer 
and CD shape on pressure wave generation and flow 
contamination. The shape of the contact discontinu- 
ity may however depend on its early history, i.e. on 
the main diaphragm rupture itself. This is another 
example of causal interaction, which forces to model 
the diaphragm rupture and subsequent evolution, in 
order to at least estimate the effect on CD shape. 
This is another difficult problem which will require 
a lot of effort: the calculation done so far is only a 
first attempt at solving the problem, which helped 
us identify the areas where further improvement is 
necessary. 

Although the opening and petaling of the steel di- 
aphragm is a fundamentally 3-dimensional process, 
it is not required at this stage to include this compli- 
cation. The problem is therefore reduced to a 2D 
axi-symmetric problem, and the diaphragm open- 
ing is reduced to a case of time varying boundary 
condition. The driver and driven tubes are grid- 
ded as two separate regions: between them, the 

boundary condition is set as a reflecting wall for 
some of the grid cells, and a patching condition for 
the grid cells within the opening. The distribution 
of patched/reflecting boundary points changes with 
time; in this first attempt, the gTids are assumed 
fixed in time, i.e. the opening proceeds in steps, one 
grid point after another. If the grid had sufficent res- 
olution in the radial direction, this would be a good 
approximation to a continuous process. Some effects 
are ignored in this approximation. First, the physical 
boundary (steel) between the two regions is changing 
in time, due to the distortion of the diaphragm (the 
petaling); its dynamics should be modeled as well, 
and this will require a major effort in the future. We 
also assumed that the diaphragm is initially vertical, 
while in reality it is to some extent assuming a hemi- 


spherical shape, due to the pressure rise in the driver 
tube. Since we are only interested at this point in 
testing the numerical capabilities, this is not a con- 
cern. 

In Figure 16 we show the temperature contours 
of the flow transient, taken at three different times 
after the start of the opening of the main diaphragm. 
At 30 //seconds, the primary shock is just past the 
contracting section of the driven tube, immediately 
followed by the contact discontinuity and a hot layer 
of driver gas. The opening at that time is still small 
(about l/6 </l of the driven tube diameter). Notice 
that between the hot, sheared gas layer and the open- 
ing is a weak shock (blue contours in the figure) is 
emanating from the diaphragm opening. This is the 
result of the step in diaphragm opening: as one more 
cell of the driver section is put in contact with the 
driven tube, the sudden change in local boundary 
condition creates this weak shock. This purely nu- 
merical effect was difficult to estimate a priori , and 
is an interesting observation in itself: it forces us to 
reconsider the technique for future computations, in 
order to have a smoother, continuous opening. At a 
later time (54 //seconds), the primary shock has trav- 
elled further down the driven tube and significantly 
straightened. Oblique shocks reflected from the walls 
of the constricting section interact near the axis to 
form a strong mach disk. The shear layer has at- 
tached itself to the walls, and the driver gas behind 
the primary shock has developed a concave shape 
near the axis. At 70 //seconds, the primary shock 
is completely planar, the mach disk has shrunk, and 
a very complex flow structure follows the primary 
shock and CD. Since the stepwise diaphragm open- 
ing has produced spurious shocks, it is difficult to 
identify the real physical effects. Presumably, this 
structure is actually simpler (at least in this 2D ap- 
proximation). We can safely assume that a few main 
conclusions can be drawn from this preliminary re- 
sult: 

• the primary shock becomes planar very rapidly 
(~ 60/iseconds). 

• a complex and unsteady flow structure, domi- 
nated by a mach disk, is formed behind the CD 

• the CD itself develops a complex shape. 

The last item can be more clearly demonstrated in 
Figure 17, which shows the CD and the main shock, 
for the same time sequence. It is clear from Figure 
17-c that the CD is definitely non-planar; indeed, it is 
likely that from that moment on, the CD evolution 
will be dominated by Rayleigh-Taylor instabilities, 
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and will not recover an ideal, planar shape. 

The remainder of the simulation, i.e. the full 
propagation of this structure down the whole length 
of the driven tube, would take a considerable amount 
of CPU, and has not yet been attempted. Future 
work will focus again on the diaphragm opening and 
on a numerical technique that more accurately de- 
scribe it. It is important also to point out that this 
calculation is difficult for an additional reason: as 
the diaphragm opens, there is a strong jet of gas 
that flows through a small opening in a ideally re- 
flecting surface. This jet has a tendency to entrain 
some flow around the opening, producing an artifi- 
cial ’cavitation’ near the surface. The same numer- 
ical phenomenon can be also observed occasionally 
near the base wall of a wedge in a hypersonic flow, 
for example. This numerical instability can be re- 
moved by enhancing the numerical diffusion in that 
region. This was done so in our case. This will also 
affect the diffusion of the driver gas into the driven 
tube, and therefore the extent of mixing and shape of 
the contact discontinuity. This transient flow must 
therefore be recomputed with special care. 

Another example of causality relation concerns 
the influence of the combustion process in the driver 
and the flow in the driven tube. Pressure fluctua- 
tions are experimentally observed in the driver tube, 
which seems to indicate that the combustion is not 
a uniform process, and that some pressure waves are 
bouncing back and forth in the driver tube. These 
waves can travel down the driven tube after opening, 
and influence the pressure field in the stagnation re- 
gion. This effect is demonstrated in Figure 18, where 
the pressure histories at two locations are compared 
with the computations: Figure 18- a shows the pres- 
sure history a few centimeters upstream of the main 
diaphragm, in the driver tube. Figure 18-b shows 
similar profiles near the end wall of the driven tube. 
Pressure fluctuations with a sine shape have been 
superimposed in the driver tube at the moment of 
diaphragm rupture; the amplitude of this fluctuation 
was chosen to match the experimental observation in 
the driver tube. These quasi- ID computations show 
a very good agreement with the experimental traces 
when the fluctuations (of the right phase) are super- 
imposed. The agreement is notably much better than 
without these fluctuations (Figure 18-b). It is clear 
that the combustion process in the driver tube has a 
strong influence on the flow conditions in the stagna- 
tion region, and should also be modeled and better 
understood. These calculations were performed by 
G. Wilson at NASA-Ames, using a quasi-ID code: 


additional details on the modeling will be presented 
in the future [7]. 

These results are another example of the useful- 
ness of dimensional reduction: it also shows that with 
some clever modeling, one can causally reduce the 
problem and retain an accurate description of the 
system. This should also help in understanding the 
conditions for better reproduction of test runs, and 
for more uniform flow conditions. This is another 
difficult task, which involves the modeling of energy 
deposition, initiation of deflagration and flame prop- 
agation. An effort in that direction is planned. 

V: Conclusions and Future Plans 

From these preliminary results, we can draw the 
following conclusions: 

• The flow conditions are very severe and put 
enormous strain on the accuracy and stability 
of the current numerical techniques. Further re- 
search into the improvement of the numerical 
techniques is desirable. 

• The validation of the CFD capabilities will re- 
quire some difficult measurements, including 
high resolution video recording of the tran- 
sient processes, such as secondary and main di- 
aphragm rupture. 

• The numerical modeling of the facility, even 
with simplistic assumptions (quasi- ID), can 
greatly benefit the design of experiments, diag- 
nostic procedures, new test conditions, and un- 
derstanding of the tunnel performance. 

• The computation of separate regions of the facil- 
ity treated as independent is an approximation 
at best. The influence of remote and past events 
on the overall flow structures is not negligible, 
and must be estimated. The separation of the 
facility flow into separate regions can be used for 
the estimation of various physical phenomena. 

The most immediate challenge still concerns the 
opening of the main diaphragm. The preliminary re- 
sults shown here suffer from an inadequate treatment 
of the unsteady boundary condition, and must be 
recomputed more accurately. Other challenges that 
await us concern the importance of various physical 
effects, such as shock-BL interactions at the end of 
the driven tube, and wall chemistry and ablation and 
resulting flow contamination. These three topics will 
be addressed in the future. 
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Figure 1: 


Schematic (not to scale) of Ames 16” Shock Tunnel facility. All dimensions are in cm. 



Figure 2: Computation of stiff chemistry. The open symbols are obtained when the implicit chemistry 
solver is sub-iterated with an increasingly large sub-step. The temperature change is also estimated at each 
sub-step. The integration proceeds until At = 5 10~ 7 seconds, which is an assumed global time step. Deviations 
from the ’exact’ solution at times greater than 10" 7 seconds are due to an error in the estimate of temperature. 
As the first integration is terminated, a more accurate estimate of the temperature (including variations of the 
specific heat) is made, and the solution quickly converges to the correct one. The insert shows the effect of 
rescaling of the chemical time scale: here the integration proceeds only to 10% of the global time step. The 
integration during the second global step also quickly converges to the exact solution. 
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Figure 3: Grid subscaling performance. A constant fixed grid (1000 points) is compared with the results 
from a high density grid (800 points) moving along with the shock. This final profile, at the end of teh driven 
tube, shows the increase in resolution of the contact discontinuity. The high density ( subscaled ) grid is a subset 
of the ’background’ grid (250 points), and moves by steps equivalent to one cell spacing of the background grid. 
The transfer of information between the two grids proceeds by volume averaging, and is fully conservative. 
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Figure 4: Driven/Nozzle flow transient at 80 
/isec after shock arrival. The shock comes from 
the left, instantaneously ruptures the secondary di- 
aphragm, propagates down the nozzle and is partially 
reflected back into the driven tube. Cells are clus- 
tered near the axis, the calculation id for inviscid, 
non-reacting flow, spatially 2 n< *-order. The conical 
structure is believed to be an artifact from the nu- 
merical method. 


11 



0 



normal component 
dS k /Vul = 


AO Ax Ar 


{AO l*l)r 2 Ax At 


[AO/ 2)v 2 Ax 


dSj/Vol : 


Ar 


AOrAr 

{A0/2)r 2 Ax 


4 

Ax 


Figure 5: Schematic of finite- volume computational cell in cylindrical geometry. Surface to volume ratios 
are independent of the wedge angle. The momentum source term in the axi- symmetric formulation of the Euler 
equations is proportional to dSk , and behaves as 1/ Ar near the axis. 




Figure 7: Effect of viscosity for low nozzle pressures. The top section shows the results for a cutoff of 
gradients at a scale comparable to 0.1 mfp. The bottom section uses a cutoff at 10 mfp. The case on the left 
is for a nozzle pressure of 100 /iTorr, the case on the right is for 100 milliTorr (the experimental condition for 
the 16” shock tunnel). The cutoff allows greater stability and greater time steps in computing implicitely the 
viscous terms. The contours shown on the right hand side show that for our case of interest, the results are 
insensitive to the choice of cutoff. 
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Figure 6: Driven/Nozzle flow transient. This computation is obtained after coarsening the grid in the 
radial direction near the axis. The shock convergence on the axis (70 psec) now leads to a stable shock structure 
(’bubble’) near the axis. The formation of an intense vortex is evident for 200 /rsec. Plots at later times would 
show that the reflected shock becomes essentially planar at ~ 500/xsec, while the vortex structure persists. 




Figure 8: Nozzle flow transient (Mach contours). The curved primary shock can be identified on the right. 
Weak oblique shocks are emanating from the walls and converge to form a strong (upstream facing) mach disk 
on the axis, and reflect back. Strongly sheared flow is visible between the primary shock and the first oblique 
shocks. A contact discontinuity behind the primary shock has been numerically diffused beyond identification. 
A small vertical break in contours is due to an interpolation error by the graphics, between two different grids. 



Figure 9: Experimental shadowgraph, taken from [6]. The upstream facing shock and the attached oblique 
shocks are clearly defined in this picture. Contact discontinuities are also very clear, and show strong vortices. 
All features of Figure 8 are contained in this picture. 
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Figure 10: Interaction in test chamber (Mach contours). In this example, as the primary shock diffracts at 
the end of the cone, the remainder of the shock structure has just reached the apex of the cone. For this case, 
steady flow over the testing body is achieved within 400 /^seconds. 
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Figure 11: Sensitivity of various flow quantities to effective nozzle area. The nozzle shape is replaced by 
a modified nozzle contour from the throat on, by removing a fixed fraction of the radius (i.e. the 87% nozzle 
removes 13% of the radius), to account for the boundary layer. Examination of the plots shows that the static 
pressure is the most sensitive variable, and therefore the most useful variable to determine experimentally. 
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Figui'e 12:Compariscm of quasi- ID solutions with 2D axi-symmeiric solution. Except for the flow within 
the boundary layer, the agreement is very good for the 87% case. Most variables have a flat profile in the 
core region, except the static pressure. Again, measurement of the static pressure should yield some valuable 
information on the nozzle flow, and for code validation. 





TEMPERATURE AT STATION N3 MOLE FRACTION OF NO AT STATION N3 



Figure 13: Temperature and NO mole fraction at two stations. There is still some significant core variation 
of the temperature at N3. The NO concentration profile remains essentially frozen, with most of the gradients 
in the boundary layer, as expected. 
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Figure 14 : Computed Intensity at station N3, with and without the boundary layer. Figure 14-a (top), 
shows the intensity in absolute value, with the boundary layer included. Figure 14-b (bottom), shows the ratio 
of computed intensity with the boundary layer, versus the computed intensity without boundary layer. A similar 
plot for the exit plane would show a ratio at peak absorption of 0.2. 
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Figure 15: Comparison of quasi-ID and axi-symmetric cases (without BL). Figure 15-a (top) shows the ra- 
tio of the computed intensity for the quasi-ID solution (87%) versus the computed intensity for the 2D solution 
at station N3. Figure 15-b (bottom) shows a similar quantity for the exit plane. 













Figure 18: Comparison between experimental and computed pressure traces at two locations. Figure 18- a 
(top) shows the comparison at a few cm upstream of the main diaphragm in the driver tube. Figure 18-b (bot- 
tom) shows the traces for a point near the end of the driven tube. Significantly better agreement is obtained 
after assuming an ad-hoc pressure gradient within the driver tube, at the end of combustion. 
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DECOUPLED PREDICTIONS OF RADIATIVE HEATING IN AIR 
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Abstract 

The radiative emission along the stagnation stream- 
line and the radiative heating at the stagnation point 
of a blunt-nosed vehicle entering the Earth’s atmo- 
sphere at hypersonic speed are estimated using a par- 
ticle simulation technique with decoupled radiation. 
The fluid mechanics of the weakly ionized flow is 
computed using the direct simulation Monte Carlo 
method, DSMC. Analysis of radiation is decoupled 
from the flowfield and is estimated using NEQAIR, 
a computer code written by Park. The results are 
compared with previous DSMC computations which 
employed a simplified, coupled radiation model. The 
effects of recent advances in modeling relaxation and 
dissociative behavior with the DSMC technique are as- 
sessed in terms of radiative emission. It is found that 
the introduction of the new models decreases the pre- 
dicted total radiative heating at the stagnation point 
of the vehicle by a factor of 15. The DSMC approach 
is also compared with a continuum flow model: in each 
case prediction of radiation is decoupled from the flow- 
field. Similar sets of vibrational and chemical relax- 
ation rates are employed in these simulations. Despite 
large differences in the computed flowfield, which ex- 
hibits strong thermo-chemical nonequilibrium, the to- 
tal predicted radiative heating estimates agree within 
a factor of 2. 

Introduction 

The motivation for the present study arises from 
the requirement for accurate radiation estimates for 
hypersonic flight vehicles. These are necessary for ad- 
equate thermal protection of the spacecraft during en- 
try into the Earth’s atmosphere. A previous study for 
the Aeroassist Flight Experiment (AFE) has been re- 
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ported by Whiting and Park 1 in which flowfield data 
were obtained using the Navier-Stokes equations at the 
lower altitudes traversed by the AFE during its sweep 
through the upper atmosphere. However, at higher 
altitudes, numerical difficulties were found in obtain- 
ing solutions. These difficulties are presumed to arise 
from the failure of the Navier-Stokes equations when 
the Knudsen number of the flow is too high. 

The primary aim of the present study is to obtain 
radiation estimates using a particle simulation method 
for low-density, hypersonic flows in which strong effects 
due to thermo- chemical nonequilibrium are present 
and to compare the results with those from Ref. 1 for 
continuum flow. This is accomplished using the di- 
rect simulation Monte Carlo method (DSMC) to pre- 
dict the one-dimensional flowfield along the stagnation 
streamline. The radiative emission is then computed 
using the computer code NEQAIR 2 with the radiation 
decoupled from the flowfield solution. 

This is the first time that radiative heating has been 
estimated from DSMC computations using the ap- 
proach in which radiation is decoupled from the flow- 
field solution. A new particle simulation code has been 
developed for the study. The code contains many re- 
cent developments which have improved considerably 
the modeling of thermo-chemical relaxation with the 
DSMC technique. 

The new code is first evaluated through comparison 
against previously published particle computations of 
hypersonic, low-density, radiating flow. These previ- 
ous DSMC computations employed simplistic thermo- 
chemistry models, and the effect of the introduction of 
the improved models is assessed. Additionally, com- 
parison is made between the new code and continuum 
techniques for flow conditions corresponding to a tra- 
jectory point of the AFE vehicle. The continuum code 
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is SPRAP (Stagnation Point Radiation Program) writ- 
ten by Park 3 which solves the Navier-Stokes equations. 
Radiative emission for these continuum calculations is 
also decoupled from the flowfield solution. 

Description of the Particle Simulation 

The particle simulation code for computing one- 
dimensional shock waves in air is highly vectorized, 
and has been extended from the implementation for a 
reacting gas described by Boyd 4 to include electrons 
and ionizing reactions. An eleven species (N 2 , N, 02, 
O, NO, N 2 + , N + , 0 2 + , 0 + , N0+, E“) real air model 
is employed. The code contains many recent model- 
ing developments which have improved the ability of 
the DSMC technique to simulate thermochemical re- 
laxation phenomena. 

The code simulates rotational 5 and vibrational 6 en- 
ergy exchange using probability functions which are 
evaluated using the energy of each collision. These 
represent an improvement on previous models in which 
constant probabilities were applied. The code employs 
the Vibrationally-Favored Dissociation model of Haas 
and Boyd 7 and a corresponding recombination model. 4 
From the DSMC simulation, it is assumed that the vi- 
brational and electronic temperatures of the gas are 
equal for the purposes qf computing the radiative emis- 
sion. Therefore, it is expected that the improved phys- 
ical models employed in the DSMC code may affect 
significantly the total radiative heating predicted. 

The introduction of electrons into the code requires 
special consideration. Due to their very low mass, 
the collision rates associated with electrons are about 
two orders of magnitude higher than those associated 
with the heavy species which occur in air. This re- 
quires the use of a numerical time-step which is two 
orders of magnitude smaller than would be employed 
were the electrons absent from the flowfield. In the 
present implementation, this problem is solved by us- 
ing a time-step based on the heavy species to move the 
particles through the flowfield, and then subdividing 
the time-step into one hundred sub-steps to perform 
the collisions. Charge-neutrality is enforced through- 
out the flowfield by forcing each electron to remain 
in the same computational cell as the ion with which 
it was initially formed. While this is physically un- 
realistic, this approach should not affect the predic- 
tion of radiative emission to any great extent. Also, 


any charged particles which collide with the cold sur- 
face of the vehicle are assumed to recombine to the 
neutral species. The special considerations described 
here for simulating charged species add a significant 
numerical overhead to the calculations (although the 
DSMC code is still highly vectorized). This overhead 
is minimized by initially running the DSMC code with- 
out the charged species until a steady state is reached. 
The ionizition and charge-exchange reactions are then 
turned on, and the system is allowed to reach a new 
steady state, before sampling of flow variables begins. 

Two different sets of chemical rate data are em- 
ployed in the DSMC code. The first corresponds to 
that used by Moss et a/.. 8 This set was implemented 
without coupled vibration-dissociation: In addition, 
for each reaction, the reverse rate expressions were de- 
termined by evaluating the electronic partition func- 
tions in the equilibrium constants at a fixed tempera- 
ture. The second set of rate constants is based on those 
employed by Whiting and Park in the continuum code 
SPRAP. In their analysis, the reverse reaction rates are 
obtained by evaluating temperature-dependent curve 
fits for the equilibrium constant. The curve fits take 
the following form proposed by Park 3 : 

K t (T) = exp(Ai + A 2 ln{z) + A^z -I- A 4 z 2 + A 5 z 3 ) 

where the A{ are constants and z=10000/T. Unfortu- 
nately, this form for the equilibrium constant is not 
mathematically convenient for implementation in the 
DSMC chemistry models. 

The goal of this part of the present study is to evalu- 
ate differences in the chemical models employed in the 
continuum and particle methods. To limit the num- 
ber of factors involved in these comparisons, it is the 
aim to maintain consistency between the relaxation 
rates employed in the solution techniques. Therefore, 
a form for the equilibrium constant which takes the 
traditional Arrhenius form is fit as a function of tem- 
perature to Park’s expression. This approach has been 
very successfully applied by Boyd and Gok^en 9 to the 
ionizing reactions for N 2 . By achieving good corre- 
spondence between the two sets of chemical rate con- 
stants, very good agreement was obtained for flow so- 
lutions to hypersonic shock-waves in N 2 using parti- 
cle and continuum methods. For air, good agreement 
is generally obtained between the new DSMC expres- 
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sions and Park’s results over the temperature range of 
interest, i.e. from 10,000 to 20,000 K. 

Description of NEQ AIR 

The code employed to estimate the radiative emis- 
sion from the flowfield solution is called NEQAIR and 
was developed originally by Park.^ NEQAIR calcu- 
lates the nonequilibrium electronic excitation of the 
species in the flow and the radiation emitted by those 
species. A principal assumption made in the code is 
that the quasi-steady state (QSS) model may be used 
to determine the population of each electronic state of 
a species after the species density is found from the 
reacting-flow solution. 

The QSS model assumes that the individual rates 
of populating and depopulating each electronic state 
are much faster than the rate of change of the pop- 
ulation of the state itself. Thus, the population of a 
state is given by the difference of two large numbers 
that are about equal. This assumption is not valid in 
the flow immediately behind the shock wave where the 
populations of the excited states are low and increas- 
ing rapidly. However, as the excited states become 
populated, the species begin to emit radiation, and by 
the time the radiation reaches a significant level, the 
assumption is reasonably valid. 

As the QSS model is applied after the chemistry 
portion of the calculation is completed, the electronic 
excitation processes are assumed to be independent of 
the chemical reactions. This simplifies the calculation 
enormously. The QSS model allows a separate effective 
electronic temperature to be defined for each excited 
electronic state. These temperatures are defined by 
comparing the population of an excited state to the 
population of the ground state. At each flowfield data 
point, the QSS model implemented in NEQAIR em- 
ploys the translational and vibrational temperatures 
from the DSMC or continuum simulation. In addi- 
tion, the QSS model assigns many temperatures for 
electronic excitation (one for each excited electronic 
state of every species). 

In the estimation of radiative emission from the 
flowfield solution, NEQAIR is executed over 12 differ- 
ent spectral regions which are listed in Table 2. These 
have been determined by Whiting and Park 1 to be 
those which make significant contributions to the total 
radiative heating for hypersonic flows of air. The fol- 


lowing molecular bands have been included: N 2 + (l*)i 

N 2 (1+), N 2 (2+). N 2 (BH1), NO/3, N0 7 , 0 2 (SR). 

The ultimate test of such codes is how well the re- 
sults compare with experiment. Park 10 has compared 
calculated radiative results with shock-tube, ballistic- 
range, and Earth-entry data covering a wide range of 
flight conditions and finds good agreement, generally 
within a factor of two. 

Results and Discussion 

This section is divided into three sub-sections. In 
the first, comparison is made between the new code 
and DSMC computations reported previously by Moss 
et al.. s A number of different comparisons are made 
to examine the ability of the new code to reproduce 
exactly the results of Ref. 8 using the same physical 
model and chemical reaction rates. Specifically, the 
DSMC code employs constant rotational and vibra- 
tional collision numbers of 5 and 50, and the chemical 
reaction rates from Ref. 8 with the degree of vibration- 
dissociation coupling set to zero. The second subsec- 
tion investigates the effects of employing the improved 
physical models. For this purpose, the DSMC code 
employs variable rotational and vibrational collision 
numbers, vibrationally-favored dissociation, and the 
new chemical rate constants derived from the contin- 
uum expressions. In the third sub-section, comparison 
is made between the DSMC code using the new models 
and the continuum computations of Ref. 1. 

Comparison with Previous Results 
Using Old DSMC Models 

The new code is first assessed through comparison 
with DSMC computations presented by Moss et at. 8 
for the AFE at an altitude of 78 km (the flow condi- 
tions are given in Table 1). The free-stream tempera- 
ture and Knudsen number were 188 K and 1.2xl0~ 3 re- 
spectively. For compatability with the study reported 
in Ref. 8, constant rotational and vibrational collision 
numbers of 5 and 50 are employed. Also, the chemi- 
cal reaction rates from Ref. 8 are used, the degree of 
vibration-dissociation coupling is set to zero, and the 
shock standoff distance is specified as 0.110 m from 
the body. A total of 1,000 computational cells and 
over 100,000 simulated particles are employed in each 
of the calculations reported in the current work. 

The temperatures computed in the present study 
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(lines) are compared with data taken from Ref. 8 (sym- 
bols) in Fig. 1. It is clear that significant thermal 
nonequilibrium effects are present. In comparison to 
the profiles from Ref. 8, the present temperatures are 
slightly higher. This effect is at least partially due to 
the absence of radiative cooling in the current work. 

The mole fractions of the major species computed in 
the present study (lines) are compared with data taken 
from Ref. 8 (symbols) in Fig. 2. The present study is in 
very good agreement with the solutions from Ref. 8. In 
addition to the general form of the computed profiles 
shown in Figs. 1 and 2, the peak electron temperature 
of about 18,000 K and the maximum electron mole 
fraction of about 0.013 are in good correspondence to 
the results of Ref. 8. 

The radiative emission along the stagnation stream- 
line given by NEQAIR for the present DSMC flowfield 
solution is shown in Fig. 3. This profile compares quite 
well with the profile given in Ref. 8, although the peak 
value is somewhat larger in the present study. Integra- 
tion of the spontaneous emission along the streamline 
to the body gives a one-dimensional radiative flux of 
64.8 kW/m 2 /sr, which is about a factor of 1.6 greater 
than the result of Moss ti al. The difference may 

partly be accounted for in terms of radiative cooling 
* 

which is omitted in the present calculations by estimat- 
ing radiative emission decoupled from the flowfield. 
Considering the significant difference in the radiation 
models employed in the two studies, the agreement is 
acceptable. 

The radiative heating at the stagnation point due to 
the radiative flux given above is found using a spherical 
cap model 1 . First, the infinite slab result is obtained 
for the total radiative heating flux at the stagnation 
point by multiplying the one-dimensional value by 2 t 
for the optically thin gas part of the spectrum, and by 
r for the self-absorbing regions. The spherical cap re- 
sult is then computed as about 80% of the infinite slab 
value 1 . This procedure gives a total radiative heating 
value of about 340 kW/m 2 , which is a very high value. 
This value is even higher than the convective heating 
rate reported in Ref. 8 as 248 kW/m 2 , 

The contributions to the total radiative heating 
made by the various spectral regions which are ana- 
lyzed are listed in Table 3 under Old Model. The rea- 
sonable agreement found between the present DSMC 


results and those obtained in Ref. 8 for the radiative 
flux, indicates that the decoupled approach to radia- 
tion assumed in this study is an appropriate solution 
method for these flow conditions. In addition, analysis 
of the rate of loss of energy from the flowfield showed 
that only about 0.2% was due to radiative emission. 
This confirms that simulation of radiation coupled to 
the fluid mechanics is unnecessary for these fiowfields. 

Comparison of Old and New DSMC models 

Having established that the new DSMC code pro- 
vides solutions which are similar to previous studies, 
when using similar models, it is appropriate to assess 
the effect on the radiative emission by the introduc- 
tion of the improved physical models, Specifically, 
variable rotational and vibrational collision numbers, 
vibrationally-favored dissociation, and the new chemi- 
cal rate constants derived from the continuum expres- 
sions are implemented. The translational and vibra- 
tional temperatures computed with the DSMC code 
using the old and new models are shown in Fig. 4 for 
the same conditions considered in the previous subsec- 
tion. The relaxation zone is much larger with the im- 
proved models. Note that the vibrational mode using 
the new model does not equilibrate with the transla- 
tional mode until the body surface is reached. This is 
due mainly to the use of the variable vibrational col- 
lision number. This quantity is a strong function of 
temperature, and only reaches a maximum of about 
50 at the peak translational temperature. For most 
of the stagnation streamline, the vibrational collision 
number is greater than 100 which reduces significantly 
the rate of equilibration of the vibrational mode. 

The species which radiate most energy in the ultra- 
violet region of the spectrum (0.1 1-0.18^ m) are atomic 
nitrogen and oxygen. The variation in the mole frac- 
tions for these species is shown in Fig. 5. With the 
new, improved models, the rise in N and 0 due to 
dissociation is retarded significantly due to the use of 
the vibrationally-favored dissociation models and the 
lower vibrational temperatures. The strongest radi- 
ator above 0.2 fim is N 2 + and its variation in mole 
fraction computed with the old and new DSMC mod- 
els is compared in Fig. 6. The rise in N? + is faster 
with the old models although behind the shock the 
two solutions show general agreement. 

The profiles of radiative emission are compared in 
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Fig. 7. The peak emission obtained with the new 
models is significantly lower than that computed with 
the previous models due to the lower electronic tem- 
perature, which is set equal to the vibrational value 
in NEQAIR. The computed three-dimensional inte- 
grated radiative heating at the stagnation point is only 
about 23.4 kW/m 2 for the new model which is a fac- 
tor of about 15 smaller than that predicted by the 
old model. The contributions to the total radiative 
heating made by the various spectral regions which 
are analyzed using NEQAIR are listed in Table 3 un- 
der New Model. The very large difference obtained in 
the radiative heating estimates using the old and new 
thermochemistry models demonstrates the sensitivity 
of the calculations to the physical models. 

Comparison with Continuum Results 

The DSMC code using the new models is applied to 
the flow conditions examined by Whiting and Park 1 , 
which are given in Table 1. These are slightly different 
from those investigated in Ref. 8. The most important 
difference involves the shock standoff distance which is 
significantly larger in Ref. 1 (0.188 m) than in Ref. 8 
(0.110 m). The freestream temperature is 188 K and 
the Knudsen number is 1.7xl0“ 3 . The DSMC compu- 
tations again employed the energy-dependent proba- 
bilities of rotational and vibrational energy exchange, 
and the vibrationally-favored dissociation model dis- 
cussed by Boyd. 4 The vibrational and chemical rate 
constants are made consistent with those employed in 
the continuum simulation as discussed earlier. 

The continuum results taken from Ref. 1 were com- 
puted using SPRAP which was developed by Park. 3 
This code computes the viscous, one-dimensional, con- 
tinuum flow behind an infinitely thin normal shock 
wave using the approach described in Ref. 3. The so- 
lution for the one-dimensional, uniform area duct may 
then be transformed to represent the flow along the 
stagnation streamline of a spherical body. The initial 
flow conditions immediately behind the shock are de- 
termined from the Rankine-Hugoniot shock-jump rela- 
tionships, assuming that the vibrational temperature 
remains at the free-stream value. A viscous shock layer 
treatment is applied for computation of the boundary 
layer which forms next to the cold wall of the vehicle. 
It is assumed that the flow in the boundary layer is 
chemically frozen with a ratio of specific heats equal 


to 11/9. Further, the wall is assumed to be noncat- 
alytic for the dissociation reactions and fully catalytic 
for the ionization reactions. In this approach, a two- 
temperature dissociation model is used in generating 
the reacting flow behind the normal shock wave. This 
model equates the molecular rotational temperature to 
the kinetic temperature of the atoms and molecules, 
and also equates the electron kinetic and electronic 
temperatures to the molecular vibrational tempera- 
ture. In the two-temperature model, all molecules 
have the same vibrational temperature, the degree of 
ionization is small, and the chemical reactions occur 
in the ground states of the atoms and molecules. 

The computed temperature profiles are compared 
in Fig. 8. The two-temperature continuum approach 
provides a translation-rotation value, and a vibration- 
electron-electronic value. The particle DSMC ap- 
proach provides separate values for the translational, 
rotational, and vibrational energy modes and equates 
the electronic temperature to the vibrational temper- 
ature to calculate the radiation. 

The temperature profiles exhibit many differences. 
The DSMC solution shows considerable shock thick- 
ness and structure which is omitted in the continuum 
calculation. In the DSMC computation, there is a sig- 
nificant region of the shock where the translational 
and rotational modes are not equilibrated, thereby 
casting doubt on the validity of the continuum two- 
temperature approach at least at these low densities. 

The rise in vibrational temperature predicted by 
DSMC is slower than that computed in the contin- 
uum simulation. Further, the DSMC method pre- 
dicts a higher degree of nonequilibrium between the 
vibrational mode and the translational and rotational 
modes all along the stagnation streamline. The rea- 
sons for these differences in the two solutions are not 
clear at this point. The factors involved include dif- 
ferences in simulating viscosity, thermochemistry, and 
transforming a one-dimensional calculation into a stag- 
nation streamline flow. These factors need to be in- 
vestigated more thoroughly. 

The mole fractions of atomic nitrogen and oxygen 
computed with the continuum method and the DSMC 
code for these conditions are compared in Fig. 9. 
Generally, the agreement is quite good except within 
the shock wave itself which is omitted in the con- 


6 AIAA 92-2971 


tinuum code. In the continuum technique, it is as- 
sumed that no chemical reactions occur upstream of 
the shock standoff location. The DSMC computation 
indicates that significant chemical activity does take 
place through the shock. This is due to the finite 
thickness of the shock-wave and to diffusion effects. 
Thus, the rise in atomic mole fractions computed by 
DSMC preceeds the continuum results by a significant 
distance. The profiles computed by the two methods 
then follow similar paths until the last few centimeters 
next to the body, where the continuum code assumes 
the flow to be chemically frozen. The variation in mole 
fraction of N 2 + computed with the two techniques is 
shown in Fig. 10. The rise in N? + predicted by DSMC 
occurs in the thick shock-front and therefore leads the 
continuum solution. Behind the shock, the DSMC re- 
sults provide a higher concentration of N 2 + than the 
continuum data. 

The spontaneous radiative emission profiles ob- 
tained from the continuum and DSMC solutions are 
compared in Fig. 11. The peak emission obtained 
with the new DSMC models is significantly lower than 
that computed with the continuum method, due to the 
smaller electronic (vibrational) temperature. 

The integrated radiative heating at the stagnation 
point computed with the DSMC calculation is about 
47 kW/m 2 while the continuum solution gives a value 
of 75.5 kW/m 2 . The contributions to the total radia- 
tive heating made by the various spectral regions are 
listed in Table 4. 

Despite the large differences observed in the flow- 
field solutions, it is interesting to note that the con- 
tinuum and DSMC total radiative heating estimates 
lie within a factor of 2 of each other. It is also worth 
noting that by increasing the shock standoff distance 
from 0.110 m (from Ref. 8) to 0.188 m (from Ref. 1) 
the total radiative heating predicted by the new DSMC 
thermochemistry models is doubled, indicating nearly 
linear scaling between these two quantities. 

For the sake of completeness, the flow conditions 
investigated in Ref. 1 were also computed with the 
old DSMC thermochemistry models. The difference 
in flow quantities computed with the old and new 
models was similar to those shown in Figs. 4 through 
6. The flow quantities were again interpretted in 
terms of total radiative heating at the stagnation point 


using NEQAIR. The solution obtained with the old 
DSMC thermochemistry models gave a radiative heat- 
ing which was ten times higher than that predicted 
with the new models. 

Concluding Remarks 

A new approach has been evaluated for predict- 
ing the radiative heating of a blunt-body entering the 
Earth’s atmosphere. In this approach, the fluid me- 
chanics of the flow along the stagnation streamline was 
predicted using a particle method (the direct simula- 
tion Monte Carlo method, DSMC), which is decoupled 
from the radiative emission. 

Comparison was made with a previous DSMC cal- 
culation which employed outdated thermochemistry 
models together with a simplified, coupled radiation 
model. For the purpose of this comparison, the present 
DSMC computation also employed the old DSMC 
thermochemistry models. These two very different ap- 
proaches gave agreement to within a factor of 2 for the 
total radiative heating at the stagnation point. This 
is viewed as acceptable considering the differences be- 
tween the radiation models employed. The decoupled 
approach gave the higher value which was attributed 
in part to the absence of radiative cooling. 

For the same flow conditions, the use of new DSMC 
thermochemistry models led to a decrease by a factor 
of 15 in the total radiative heating at the stagnation 
point. This drastic reduction in heating illustrates the 
sensitivity of the computed data to the physical models 
employed in the simulations. It is proposed that the 
new thermochemistry models, which have been suc- 
cessfully evaluated in previous studies, should provide 
more realistic simulation results. This decrease in ra- 
diative heating indicates that it is not necessary to cou- 
ple radiation to the fluid mechanics under these flow 
conditions. To attain greater confidence in the numer- 
ical simulations, it is clear that experimental data is 
required for the validation of these physical phenom- 
ena. 

A further comparison of the new DSMC chemistry 
models with a continuum calculation gave agreement 
within a factor of 2 for the total radiative heating, 
despite considerable differences in the flowfield solu- 
tions. The main conclusion of this study is that there 
remains a large degree of uncertainty in the applica- 
tion of state-of-the-art numerical methods for the pre- 
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diction of radiative heating to hypersonic vehicles fly- 
ing in low-density regions of the Earth’s atmosphere. 
Considering the importance of such heating to many 
future aerospace projects, there is a need for continued 
numerical and experimental investigation in this area. 
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Table 1. Freestream conditions. 


Source Altitude 

Uoc (m/s) 

Poo (kg/m 3 ) 

Ref. 8 78 km 

9110 

2.75xl0~ 5 

Ref. 1 80 km 

9756 

2.05xl0 -5 

Table 2. Spectral regions analyzed with NEQAIR. 

Region (pan) 

Description 

Absorption 

0.1130-0.1140 

N atomic lines 

Yes 

0.1160-0.1170 

N atomic lines 

Yes 

0.1170-0.1180 

N atomic lines 

Yes 

0.1195-0.1205 

N atomic lines 

Yes 

0.1240-0.1250 

N atomic lines 

Yes 

0.1294-0.1314 N and O atomic lines Yes 

0.1315-0.1325 

N atomic lines 

Yes 

0.1323-0.1333 

N atomic lines 

Yes 

0.1405-0.1415 

N atomic lines 

Yes 

0.1490-0.1500 

N atomic lines 

Yes 

0.1740-0.1750 

N atomic lines 

Yes 

0.2000-2.0000 

Molecular lines 

No 

Table 3. Total radiative heating (in 

kW/m 2 ) using 

DSMC for the flow conditions from Ref. 8. 

Region (/im) 

Old Model 

New Model 

0.1130-0.1140 

3.11 

0.08 

0.1160-0.1170 

4.19 

0.16 

0.1170-0.1180 

2.30 

0.08 

0.1195-0.1205 

6.07 

0.21 

0.1240-0.1250 

3.65 

0.14 

0.1294-0.1314 

8.23 

0.24 

0.1315-0.1325 

3.16 

0.23 

0.1323-0.1333 

2.37 

0.07 

0.1405-0.1415 

3.13 

0.28 

0.1490-0.1500 

11.94 

0.72 

0.1740-0.1750 

19.30 

2.21 

Sub-Total 

67.5 

4.42 

0.2000-2.0000 

272.2 

19.0 

Total 

339.7 

23.4 


Table 4. Continuum and DSMC contributions (in 
kW/m 2 ) of the various spectral regions to the total 
radiative heating at the stagnation point for the flow 
conditions from Ref. 1 at an altitude of 80 km. 


Region (/im) 

Continuum 

New Model 

0.1130-0.1140 

0.44 

0.20 

0.1160-0.1170 

0.76 

0.50 

0.1170-0.1180 

0.32 

0.27 

0.1195-0.1205 

0.84 

0.46 

0.1240-0.1250 

0.92 

0.40 

0.1294-0.1314 

0.76 

1.13 

0.1315-0.1325 

0.76 

0.61 

0.1323-0.1333 

0.52 

0.19 

0.1405-0.1415 

0.84 

0.70 

0.1490-0.1500 

2.84 

1.65 

0.1740-0.1750 

5.36 

4.73 

Sub-Total 

14.4 

10.8 

0.2000-2.0000 

61.1 

36.2 

Total 

75.5 

47.0 
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Fig. 1. Temperature profiles along the stagnation 
streamline: comparison of present calculatons and 
Ref. 8 using the old DSMC thermochemistry models. 



Fig. 2. Profiles of atomic species mole fractions: com- 
parison of present calculatons and Ref. 8 using the old 
DSMC thermochemistry models. 



Fig. 3. Profiles of radiative emission along the stagna- 
tion streamline: comparison of present calculatons and 
Ref. 8 using the old DSMC thermochemistry models. 



Fig. 4. Temperature profiles along the stagnation 
streamline: comparison of the old and new DSMC 
thermochemistry models. 



Fig. 5. Profiles of atomic species mole fractions along 
the stagnation streamline: comparison of the old and 
new DSMC thermochemistry models. 
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Fig. 6. Profiles of N 2 + mole fraction along the stagna- 
tion streamline: comparison of the old and new DSMC 
thermochemistry models. 
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Fig. 7. Profiles of radiative emission along the stagna- 
tion streamline: comparison of the old and new DSMC 
- thermochemistry models. 
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Fig. 8. Temperature profiles along the stagnation 
streamline: comparison of continuum and DSMC cal- 
culations using the new thermochemistry models. 
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Fig. 9. Profiles of atomic species mole fractions: com- 
parison of continuum and DSMC calculations using 
the new thermochemistry models. 
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Fig. 10. Profiles of N 2 + mole fractions: comparison 
of continuum and DSMC calculations using the new 
thermochemistry models. 



o 10 20 30 


Distance from body (cm) 

Fig. 11. Profiles of radiative emission: comparison 
of continuum and DSMC calculations using the new 
thermochemistry models. 
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Abstract 

The purpose of this paper is to present a numer- 
ical methodology for simulating the time-dependent 
reacting flow inside the entire length of several types 
of high enthalpy pulse facilities. The numerical ap- 
proach uses a finite- volume TVD scheme for the quasi- 
one-dimensional Euler equations coupled with finite- 
rate chemistry on a moving mesh. A Riemann solver 
is incorporated for tracking contact discontinuities. 
Several numerical difficulties encountered in simulat- 
ing pulse facility flows are discussed and the numeri- 
cal techniques chosen to overcome these problems are 
presented. Comparisons of experimental and compu- 
tational pressure traces are made for three different 
experimental facilities; the NASA Ames electric arc 
driven shock tube facility (from cold driver shots), the 
NASA Ames 16 inch combustion driven shock tun- 
nel, and the HYPULSE expansion tube. The com- 
parisons help validate the numerical methodology and 
add to the understanding of the flow inside these facil- 
ities. Differences between the quasi-one-dimensional 
solutions and the experiments are used to suggest im- 
provements that multi-dimensional simulations might 
provide. 

Introduction 

The high enthalpy pulse facilities considered in this 
work are experimental apparatus designed to create 
hypersonic flow. The term “pulse” is used because 
these facilities are capable of creating the desired flow 
conditions for only short periods of time (on the order 
of milliseconds). The present work investigates shock 
tubes, shock tunnels (without pistons), and expansion 
tubes which represent only a few of the many differ- 
ent types of pulse facilities. At the core of each of 
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these facilities is a shock tube, which makes the cur- 
rent methodology applicable to several devices. 

A primary motivation for simulating the flow in- 
side a high enthalpy pulse facility is to help predict 
important operating characteristics such as the stag- 
nation properties and the chemical composition of the 
test flow which are often difficult to measure. Know- 
ing these quantities is important for interpreting the 
experimental data obtained from these facilities. For 
example, the shape/angle of shock waves or the ig- 
nition properties in a combustor can be affected by 
test gas that has become contaminated or dissociated 
during the operating sequence of a facility. Numerical 
simulations that are able to quantify such effects can 
make the experimental data more valuable. They can 
also be used to investigate ways to increase a facility’s 
performance or operating envelope. 

Computing the unsteady flow inside a high enthalpy 
pulse facility is challenging because the flow can be sig- 
nificantly influenced by strong shock waves, gas inter- 
faces, high temperature effects, chemical kinetics, and 
viscous interactions. In addition to these complexi- 
ties, there is often a large disparity of length scales. 
For example, a facility may be tens of meters in length 
while the length scales of some of the physical phenom- 
ena such as the chemical kinetics may be less than a 
millimeter. This requires findings way to resolve the 
small scale physics while still economically simulating 
the large scale flow. 

Because of the large length scales of the facilities 
being considered in this work (in one case over 50 m), 
an unsteady quasi-one-dimensional compressible flow 
formulation has been adopted. The lower computa- 
tional requirements of the one-dimensional approach 
compared to multi-dimensional simulations allows the 
computations to be advanced over long distances in a 
practical amount of computer time and also allows for 
the use of a relatively complicated gas model. 

Several recent works have also used one-dimensional 
simulations to look at high enthalpy pulse facilities. 
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Groth ei al. 1 used a sophisticated gas model to com- 
pute the complete unsteady performance of an impulse 
tunnel. Their work clearly shows the usefulness of one- 
dimensional simulations of a complete facility. How- 
ever, that work differs from the present research be- 
cause instead of having gas interfaces, the facility that 
they simulated uses a piston to separate the driver gas 
and test gas. As a result, several of the primary phys- 
ical phenomenon and numerical concerns are differ- 
ent. Another work by Maus et al 2 contains quasi-one- 
dimensional simulations of a free-piston shock tunnel 
using a Lagrangian scheme. 

The numerical approach taken in this work has 
been to solve the unsteady quasi-one-dimensional Eu- 
ler equations coupled to a detailed finite-rate chem- 
istry model. An upwind TVD scheme is used to com- 
pute the bulk of the flowfield and a Riemann solver 
is included for tracking gas interfaces. This hybrid 
scheme is applied on a moving mesh. One goal of 
this work is simply to investigate whether this kind 
of approach can consistently predict some of the de- 
tailed flow features within large facilities that operate 
by rupturing diaphragms. There are concerns that the 
deformation and instabilities of contact discontinuities 
are so severe and unpredictable that the current ap- 
proach with its simplifications and ideal treatment of 
the gas interfaces may not be helpful. With this over- 
all objective in mind, some of the finer points of the 
physical modeling have been kept relatively simple and 
show what is possible rather than what is best. For in- 
stance, viscosity and heat conduction models are not 
as sophisticated as they might be and the provision 
for vibrational nonequilibrium is included in the fluid 
equations but is not utilized. 

This work begins by presenting and discussing sev- 
eral numerical difficulties, which are considered to be 
common to simulations of high enthalpy pulse facili- 
ties. This is followed by a description of the mathe- 
matical and numerical formulation. Then, numerical 
computations are compared to experimental data from 
a shock tube, a shock tunnel, and an expansion tube. 
In addition to validating the numerical approach, these 
simulations demonstrate the physical understanding 
that can be gained with quasi-one-dimensional model- 
ing. Specifically, the computations are used to predict 
the following phenomena: the influence of initial con- 
ditions on the interaction of the reflected shock with 
the contact discontinuity (tailoring behavior); the in- 
fluence of pressure non-uniformities in the driver sec- 


tion (a situation often seen in combustion drivers) on 
the reflected shock region; and the influence of a finite 
diaphragm opening time on the chemical composition 
of expansion tube test conditions. 

Numerical Issues 

Two numerical issues encountered in this work are 
discussed in this section so that the reasons for adopt- 
ing the present numerical scheme can be understood. 
One issue is the performance of numerical schemes on a 
flow started by instantaneously removing a diaphragm 
which separates two regions with a very large pressure 
difference, and the other is the smearing of gas inter- 
faces by numerical diffusion. 

All the pulse facilities considered in this work begin 
by bursting a diaphragm which separates a high pres- 
sure driver section from a low pressure driven section. 
One-dimensional simulations of the unsteady flow cre- 
ated by this type of initial condition are not uncom- 
mon and plots of the resulting shock, contact discon- 
tinuity, and expansion fan are often presented in the 
literature when the merits of a numerical scheme are 
being investigated. 3,4,5 These investigations, however, 
are usually for relatively low pressure ratios across the 
diaphragm (ratios less than 10). Initial pressure ra- 
tios in a hypersonic pulse facility may be over 10000 
in which case it has been found that special care must 
be taken. 

Figure 1 shows temperature profiles from simula- 
tions of a perfect gas in a 2 meter shock tube which 
began with driver conditions of 1 atm and 300K and 
driven conditions of .1 atm and 300K (initial pressure 
ratio 10). The computations were made using Harten 
and Yee’s upwind TVD scheme 4 at CFL number 0.65 
on several different equally-spaced grids. All of the so- 
lutions are similar with the major difference being the 
degree of resolution of the contact discontinuity (and 
to a smaller extent the sharpness of the shock). The 
similarity of the solutions and their convergence to- 
ward a single solution gives confidence in the accuracy 
of the computations. In contrast, Figure 2 depicts tem- 
perature profiles resulting from the same initial con- 
ditions except that a driver pressure of 1000 atm was 
used (initial pressure ratio 10000). Under these con- 
ditions, the solution changes significantly with each 
refinement of the mesh and even a 3200 point grid is 
unable to support an accurate solution. In this case, 
not only does the resolution of the features change with 
the grid, but the shock propagation speed changes as 
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well. Investigations with several first and second-order 
schemes including a Riemann solver resulted in similar 
results with the conclusion that higher grid resolution 
is required at the diaphragm as the initial pressure 
ratio is increased. 

This conclusion is supported by using a highly clus- 
tered grid at the diaphragm for just the initial part of 
the solution. To avoid the small time step and high 
cost that a fine spacing imposes, the grid spacing is 
allowed to expand as the solution progresses until an 
equally-spaced grid is attained. The temperature pro- 
files with an initial spacing of 2.5 /im at the diaphragm 
location (a 400,000 point equally-spaced grid would be 
required to get this spacing) are shown in Figure 3 
for grids of 200, 400, and 3200 points along with the 
equally-spaced 3200 point solution from Figure 2. The 
plot shows that with the initial clustering the propa- 
gation speed of the various features becomes nearly 
the same for all the grids and that a very accurate 
solution is reached more quickly than on the equally- 
spaced grid (it is not shown but an 800 point solution 
is nearly the same as the 3200 point solution). 

The other numerical issue mentioned above is the 
smearing of gas interfaces by numerical diffusion. The 
smearing increases as the grid becomes coarser and 
also increases as the distance the interface travels be- 
comes greater. The latter property is a concern when 
the length of the facility being simulated becomes 
large. When an interface is poorly resolved, the result- 
ing artificial diffusion can lead to nonphysical spikes or 
valleys near the interface (see Figure 2) and if species 
concentrations change across the interface, nonphysi- 
cal mass transfer can occur. Such behavior can pose 
additional problems if chemical kinetics are being in- 
cluded in the simulation because nonphysical chemical 
reactions may occur. 

Although some numerical algorithms produce less 
interface smearing than others, all are less than satis- 
factory without the presence of a fine grid in the region 
of the interface. However, this increases the mesh size 
and reduces the time step. Another way to avoid the 
problem is to track the contact discontinuity and elimi- 
nate the interface smearing completely. This approach 
is particularly well suited to one-dimensional problems 
because the interface cannot become distorted and so 
the tracking algorithm can remain relatively simple. 
Given these observations, a tracking scheme was de- 
veloped for this work. The scheme utilizes the same 
TYD method for all of the flowfield except at gas in- 


terfaces where a Riemann solver is used. The details 
of the scheme are given later. 

Results using the tracking scheme on the case with 
a pressure ratio of 10000 are presented in Figure 4 on 
grids of 200 and 400 points. Also included in the figure 
is the well resolved 3200 point solution from Figure 3. 
The solutions obtained with the tracking scheme also 
utilize initial clustering at the diaphragm since the 
tracking algorithm does not remove the requirement 
for high resolution at this location when there are high 
pressure ratios. As before, the minimum spacing is 
increased as the solution proceeds. The effectiveness 
of tracking the contact discontinuity is clearly seen. 
There is no smearing of the contact discontinuity and 
both the 200 and 400 point solutions closely match the 
3200 point solution without tracking. The points are 
plotted for the 200 point solution showing that there 
are no points within the gas interface. The scheme 
with front tracking required a reduced time step for 
the first 50 time steps but otherwise was run at the 
same CFL number. 

Not all of the improvement from the interface track- 
ing seen in Figure 4 is due to the removal of numerical 
diffusion across the contact discontinuity. Some of the 
improvement is simply due to increased grid resolution 
at the interface. By moving the mesh with the inter- 
face, the initial clustering at the diaphragm follows 
the interface. Therefore, a side benefit to tracking the 
contact discontinuity is that it provides an easy way 
to concentrate grid points where they are needed. 

Governing Equations 


The fluid dynamic equations for the quasi-one- 
dimensional Euler equations with chemical and vibra- 
tional nonequilibrium and approximations for viscous 
and heat-conduction losses are written in the following 
form: 
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The variable S represents the cross sectional area of 
the duct and is a known function of x. The other terms 
are the density of each species pi with p = Ya=i A‘> 
the pressure p, the velocity components u and t/, the 
vibrational energy per unit volume E Vi and the total 
energy per unit volume E. The source terms repre- 
sent the production or destruction of species i through 
chemical reactions, and the term w v is the source term 
for vibrational energy. The source term A is due to 
area changes and F wa u and Q wa ii represent the losses 
in momentum and energy due to viscosity and heat 
conduction per unit volume, respectively. 

The pressure is given by the equation of state for a 
thermally perfect gas mixture, 


n 3 



where R is the universal gas constant and Mi is the 
molecular weight of species i. The total energy is writ- 
ten as 


ns ^ ru 
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(4) 

In Equation (4), the first term is the sum of transla- 
tional and rotational energies per unit volume for each 
species ( c v is the translational and rotational specific 
heat). The second term E v + t \ contains the sum of 
the vibrational energies for all the species containing 
vibrational modes plus the sum of the energy in the 
excited electronic modes. The vibrational and elec- 
tronic modes are assumed to be in thermal equilibrium 
at a vibrational-electronic temperature T v . The vibra- 
tional energies are computed assuming harmonic oscil- 
lators with corrections for anharmonic effects applied 
to species that may be present in large mass fractions 
and have significant anharmonic terms (iV 2 , 0 2l H 2i 
and H 2 0). Electronic modes are included for only 0 2 , 
0 ,and N. The last two terms in Equation (4) are the 
kinetic energy and the chemical energy per unit vol- 
ume. The thermodynamic data are from the JANAF 
tables. 6 

The vibrational source term w v is in the form of 
Landau-Teller relaxation. In this work, very fast re- 
laxation rates are specified so that the gas is essentially 
in thermal equilibrium (i.e. T v = T) at all times. The 
source term due to area changes A is written 


A = 


pdS 
Sdx ’ 


The source term modeling the boundary layer fric- 
tional loss F wa it is written in the following form 

4/ pu\u\ 

Fwall — ~ £)"~2 

where / is a friction coefficient and D is the tube di- 
ameter. The rate of heat transfer to the gas from the 
wall Q wa u is given by 

Qwall — ~^{F wa u T) 

where T wa u is the wall temperature and a is a heat- 
transfer coefficient. The expressions for F wa u and 
Qwau are very simple and the coefficients in this work 
are adjusted by making comparisons to experiment for 
one case and then using those coefficients for all other 
cases. A more complex treatment of these terms is 
found in Groth ei a/.. 1 In whatever way these terms 
are defined, they are still engineering approximations 
for multi-dimensional phenomena and probably need 
to be “tuned” for each particular facility that is in- 
vestigated. The application of these kinds of simple 
models to flows that have large boundary layers in re- 
lation to the tube diameter may not work at all (i.e. 
a satisfactory model for the low density, high velocity 
flow in the acceleration section of an expansion tube 
has not been found for this work). 

Chemistry Model 

The chemistry mechanism is based on the hydrogen- 
air combustion mechanism (including the species 
and reactions within the hydrogen peroxide and 
NO x supplements) recommended by the NASP rate 
committee. 7 The hydrogen species are included for 
the purpose of simulating facilities with combustion 
drivers. In addition, helium is included since this light 
gas is often used in pulse facilities. The complete set 
includes 14 species ( He, N 2l 0 2i # 2 , NO, OH, N0 2 , 
H NO , HOo , H 2 O t H 2 0 2 , N , O, H) and 32 reactions, 
each of which may proceed forward or backward. He- 
lium is added as an inert species with a third body 
efficiency in termolecular reactions of 1.0. An air re- 
action mechanism is obtained by setting the species 
containing hydrogen to negligible levels. 

The equilibrium constants for each of the reac- 
tions are calculated from standard Gibbs free energy 
for each species which, in turn, are obtained from 
the NASA Lewis polynomials for thermochemical data 
found in Appendix 2 of the NASP rate committee re- 
port. These polynomials are valid from 200 to 6000 K. 
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Plots extrapolated to 7000 K show that the polynomi- 
als are still well behaved out to this temperature; how- 
ever, curve fits to higher temperatures are necessary. 
Note that these curve fits are used only for computing 
the equilibrium constants for the purpose of finding 
backward rate constants. Internal energy values do 
not rely on these curves, instead they are computed in 
the manner explained in the previous section. 

Numerical Method 

The gas dynamic equations are solved by using an 
explicit finite-volume form of the Harten-Yee upwind 
TVD scheme 4 which gives second-order spatial and 
temporal accuracy. Tracking of gas interfaces and 
a moving grid are incorporated for the reasons de- 
scribed in the section on numerical difficulties. The 
combined TVD/tracking scheme is relatively straight- 
forward. The upwind TVD algorithm is used to calcu- 
late the flux across all of the interior cell sides except 
for at gas interfaces where a Riemann solver is applied 
(see Figure 5). Tracking other features such as shocks 
is not necessary because modern flux-splitting or flux- 
differencing methods capture these features quite well. 
The Riemann solver provides the local fluid velocity at 
an interface so that the cell side associated with the 
interface can be moved with it (i.e. it is tracked). As 
the interface moves, all the other points in the com- 
putational domain are moved in an organized way so 
that a distributed mesh is retained. Several minor 
changes are required to the TVD algorithm to account 
for changes in the inviscid flux due to the motion of 
the grid. Note that with the approach just outlined, 
a contact discontinuity behaves as any other interior 
part of the simulation (i.e. waves are allowed to pass 
through the interfaces). 

Because it is prescribed that there is no mass flux 
across an interface, only two quantities need to be cal- 
culated by the Riemann solver. These two quantities 
are the local fluid velocity and the local static pressure. 
The inviscid flux vector at the contact discontinuity 
F c d is written: 

F crf = (0,0,...,0,^0,pV) (5) 

where p* and u* are supplied by the Riemann solver 
and the cell side is translated at the velocity u*. These 
two quantities are generally the first two quantities 
computed by such a solver and only constitute a frac- 
tion of the computational work usually associated with 


such a solver. Therefore, there is little, if any, increase 
in cost introduced by the interface tracking. 

The procedure for the Riemann solver to obtain u* 
and p* was taken from Jacobs 8 and generalized for 
a multicomponent mixture. A MUSCL interpolation 
is used to achieve a higher order spatial accuracy. A 
modification is made in the TVD scheme for the cell 
sides adjacent to a gas interface so that the computa- 
tional stencil does not cross the interface. The effect 
of the Riemann solver, which is first-order in time, on 
temporal accuracy has not been formally established. 
No noticeable change in accuracy has been observed by 
incorporating it with the second-order upwind TVD 
algorithm. Tests of temporal accuracy were made by 
comparing solutions with and without tracking at var- 
ious CFL numbers. 

The motion of the computational mesh is dictated 
by both the motion of gas interfaces and the require- 
ment to have very fine grid initial spacing at di- 
aphragm locations. Because interfaces are tracked, 
cells initially associated with a region of gas (such as 
the driver or driven gas) remain associated with that 
region. Therefore, when the driven section is com- 
pressed by the incident shock, the cells become com- 
pressed too. This has the positive result of producing 
a fine mesh at the reservoir region of a shock tube or 
shock tunnel. However, in some cases, the grid spacing 
can become so small (and the time step so small) that 
the solution becomes very expensive. In such cases, 
the difficulty can be avoided by turning off the track- 
ing and the grid motion (reverting to the TVD scheme 
on a fixed grid for all cells) when a reasonably fine grid 
has been created. In this way, the driven section grid 
is allowed to become fine enough to satisfactorily “cap- 
ture” any gas interfaces but not so fine as to make the 
time step unpractical. 

The source term due to the area change is treated 
in an explicit manner, however, the source terms rep- 
resenting the finite-rate chemical kinetics and vibra- 
tional relaxation are often large and make the algo- 
rithm too stiff to be advanced explicitly. To avoid 
this stiffness, a fully-coupled point-implicit treatment 
of these source terms is implemented. The extra cost 
of inverting a singe block matrix at each cell is more 
than offset by the increase in the size of the time step 
that this allows. 

Shock Tube Simulation 

The first pulse facility considered is the NASA Ames 
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electric-arc driven shock tube facility. 9 * 10 This facility 
has several different configurations and has a large hy- 
pervelocity operating range using its arc-driver; how- 
ever, in this work, several shots using a cold helium 
driver with nitrogen in the driven section are simu- 
lated. These shots were made with a cylindrical driver 
of .76 m (2.5 feet) in length and 10 cm (3.93 inches) in 
diameter. The driven section was 4.02 m (13,17 feet) 
long with the same diameter as the driver. A single 
self-break diaphragm separated the driver and driven 
gases. 

Three different run conditions are simulated. All of 
the shots had a nominal 3.96 MPa (575 psi) cold he- 
lium driver. The initial nitrogen driven pressures were 
different for each shot. Driven tube pressures of 300 
torr, 100 torr, and 20 torr were used which resulted 
in incident shock speeds of approximately 1250 m/sec, 
1550 m/sec, and 1950 m/sec, respectively. The ini- 
tial conditions prescribed for the simulations were the 
same with the additional specification of a driver and 
driven tube temperature of 300 K (the experimental 
driver pressure was not known). The temperatures be- 
hind the reflected shock range between 1400-4500 K, 
so the vibrational excitation and dissociation of the 
nitrogen driven gas is important. Figure 6a-c shows 
experimental and computed sidewall pressure traces 
2.5 cm upstream of the driven section endwall. Each 
of the plots has the same scaling so that the relative 
changes caused by the different driven tube pressures 
can be easily identified. The source terms for viscos- 
ity or heat conduction were turned off to demonstrate 
what kind of results can be expected from purely in vis- 
cid one-dimensional flow. Each of the simulations used 
350 points (half in the driver and half in the driven sec- 
tion) and required approximately 5 minutes of CPU 
time on a single processor on the NASA Ames Cray 
Y-MP. 

Several of the major events are labeled in the Fig- 
ure 6. These include the passage of the incident and 
reflected shocks, the arrival of the rarefaction wave 
(which has reflected off the opposite end of the facil- 
ity), and the presence of reflected waves off the contact 
discontinuity. Agreement between the experiments 
and the computations is relatively good. Certainly, the 
general changes in the pressure traces due to reflected- 
shock/contact-discontinuity interaction (in off-tailored 
conditions) are captured for all three cases. As ex- 
pected, the experimentally measured features associ- 
ated with the the contact discontinuity (such as waves 


which have interacted with it) are smeared because 
the interface is non-planar and diffuse while the same 
features in the simulations are sharp because the in- 
terface is tracked. The position of the pressure trans- 
ducer within the boundary layers further contributes 
to the smearing of the experimentally observed fea- 
tures. Shock/boundary-layer interaction of the re- 
flected shock is large and the one-dimensional formu- 
lation cannot capture such effects. The overshoot just 
after the arrival of the reflected shock is from the 
lambda shock formed at the foot of the reflected shock 
(see Figures 6b,c). The 20 torr driven case is expected 
to be influenced most by boundary layer effects be- 
cause the shock speeds are highest and the density 
is the lowest. This is the case in which the largest 
differences are seen between the experiment and the 
computation. 

As an aid to understand the features seen in the 
pressure traces, a computed x-t diagram of density 
contours for the 20 torr case is presented in Fig- 
ure 7. Note that a vertical line of data from this plot 
yields a trace at a given location. Included in the fig- 
ure is an enlargement of the driven section endwall. 
The enlargement shows the reflected shock being par- 
tially transmitted through the contact discontinuity 
and partially reflected back toward the endwall. This 
occurs several times and each successive interaction 
further slows the interface down until it is brought 
to rest. This type of off-tailored interaction (with re- 
flected shocks) is referred to as an overdriven case. 

Shock Tunnel Simulation 

The next facility considered is NASA Ames 16-inch 
combustion-driven shock tunnel. u ’ 12 In a shock tun- 
nel the high-pressure, high-temperature region of stag- 
nated gas (created at the end of a shock tube by the 
reflection of the incident shock) is expanded through a 
nozzle to obtain supersonic/hypersonic flow. Figure 8 
contains a schematic of the facility. The driver and 
driven tubes are 21.3 m (70 feet) and 25.9 m (85 feet), 
respectively. The inside diameter of the driver tube 
is 43.18 cm (17 inches) and the inside diameter of the 
driven tube is 30.48 cm (12 inches). The name of the 
tunnel is derived from the 16-inch naval guns which 
were used to construct the driver. The nozzle is 5.9 
m (19.35 feet) in length and the Mach number at the 
nozzle exit varies between 5 and 7. 

The driver tube is operated by burning a mixture of 
hydrogen and oxygen with helium or nitrogen added as 
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a diluent. The mixture is ignited by heating thin wires 
strung down the length of the driver with a voltage 
discharge from a capacitor bank. Post-burn pressures 
are typically between 300 and 600 atm. A single self- 
break diaphragm is used which is designed to rupture 
near the peak of the pressure rise. The driven section is 
usually filled with air between 100 and 300 torr which 
provides incident shock speeds around 3 km/sec. The 
high pressure and temperature reservoir created when 
this shock reflects off the endwall lasts for 10-30 msec 
and can produce 2-10 msec of uncontaminated flow 
depending on the conditions. 

The simulations presented for this facility will con- 
centrate on the effect that the driver conditions can 
have on the reservoir. Many combustion driven pulse 
facilities have found that it is very difficult to ex- 
perimentally produce uniform burning of the mix- 
ture which then causes non-uniform post-burn driver 
conditions. 13 Sometimes these non-uniformities can be 
severe; however, a burning procedure has been devel- 
oped for the NASA Ames shock tunnel that produces 
only mild pressure oscillations (10-15% of the post- 
burn pressure) that are reproducible in amplitude and 
phase from shot to shot. 

Figure 9 shows pressure traces from a relatively low 
pressure condition in which the diaphragm was pur- 
posely kept from breaking (i.e. a closed driver shot). 
The traces correspond to the three pressure transducer 
locations in the driver depicted in Figure 8 (DR1, DR2, 
and DR3). They all show a pressure rise due to com- 
bustion over the first 20 ms followed by oscillations at 
the post-burn pressure. The traces reveal that when 
the pressure reaches a maximum at one end of the 
tube, it is at a minimum at the other end. The pres- 
sure in the center of the driver remains fairly constant. 
This type of behavior is caused by the driver gas trans- 
lating in one direction down the tube, coming to a rest, 
and then translating back again. 

The pressure traces in the closed driver experiment 
were used to find an approximate model for the flow in 
the driver. It was found that at the point in time when 
the diaphragm region (DR3) experiences a pressure 
maximum, that the driver flowfield can be modeled by 
a sinusoidal pressure and density distribution (with 
a constant temperature) superimposed on the average 
post-burn conditions as plotted in Figure 10. The fluid 
is assumed to be at rest throughout the driver at this 
point and the amplitude of the sine wave is deduced 
from the experiments. Using the approximate flowfield 


shown in Figure 10 as initial conditions for a closed 
driver simulation results in the pressure traces shown 
in Figure 11. These computed traces have good qual- 
itative agreement with the experimental traces shown 
in Figure 8. The pressure waves steepen with time 
and the middle trace has twice the frequency of the 
end traces. 

Having established some of the qualities of the 
driver, the conditions and results for an particular shot 
are presented. The case of interest used a pre-burn 
driver mixture of 

2H 2 + 1.170a + 9.36tfe 

at a pressure and temperature of 740 psi and 300 K. 
The driven section was filled with air at 220 torr. Al- 
though the pre-burn conditions are well known, the 
simulations require the post-burn pressure, tempera- 
ture, and species mass fractions as the initial condi- 
tions. Obtaining these conditions by simulating the 
time dependent ignition and combustion of the mix- 
ture is a multi-dimensional problem beyond the ca- 
pabilities of a one dimensional formulation. There- 
fore, post-burn equilibrium conditions corresponding 
to the pre-burn mixture are found using a code called 
STANJAN. 14 By this method, the ideal post-burn 
pressure and temperature for the mixture above are 
430 atm (6300 psi) and 2690 K. As a correction to the 
ideal conditions, a pressure and temperature of 90% of 
these values were used in the simulation. This partic- 
ular loss value has been found to work well over a wide 
range of facility operating conditions and accounts for 
losses due to phenomena such as heat-conduction dur- 
ing the 20 ms burning process, lateral waves, and en- 
ergy absorbed by the diaphragm xupture. 

Figure 12a, b show experimental and computed pres- 
sure traces at locations just upstream of the diaphragm 
(DR3) and next to the driven tube endwall (DN5). 
The figure includes one simulation that began with a 
uniform driver and one that began with the sinusoidal 
pressure distribution shown in Figure 10. This initial 
sinusoidal pressure distribution with a maximum oc- 
curring at the diaphragm location when the diaphragm 
ruptures was deduced from the experimental traces 
which show that this is almost always the case (it is 
logical that the diaphragm breaks while experiencing 
a pressure peak). The time basis for the simulations 
and the experiment were aligned by placing the start 
of the simulation at the pressure drop corresponding 
to the diaphragm rupture at DR3 (at approx. 29 ms). 
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It can be seen that the initial conditions that include 
the driver pressure gradients give a much better over- 
all agreement with the experiment and reveal how the 
driver pressure oscillations effect the pressure in the 
reservoir. Such a finding would be difficult without a 
simulation of the whole facility. 

The sharp drop and then rise in pressure seen at 
DE3 immediately after the diaphragm break is due to 
the rapid area change between the driver and driven 
tube (the area changes by almost a factor of 2). The 
drop exists in both the experimental and simulated 
traces but is much larger in the computations. The 
large differences between the experiment and simula- 
tions may exist in the diaphragm region because of 
two dimensional and viscous effects caused by the area 
change. It was found that having rapid area changes 
increased the required number of grid points for a 
smooth solution. These simulations used 600 points 
and required approximately 40 minutes of CPU time 
on the NASA Ames Cray Y-MP. Solutions with 1200 
points gave nearly the same results. 

The viscous coefficient was set so that the simu- 
lations would duplicate the attenuation of the inci- 
dent shock as measured by pressure transducers lo- 
cated down the length of the driven tube. It was found 
that a single value of this coefficient was adequate for 
a variety of facility operating conditions. The coeffi- 
cient for heat loss was set so that closed driver simu- 
lations would predict a pressure loss over time similar 
to that observed in experiments (see Figure 9). Only 
one value of the heat-conduction coefficient was used 
for a variety of conditions. With the heat-conduction 
coefficient set in this way, the simulations showed that 
heat-conduction was not a relatively important aspect 
of the modeling for this particular facility. 

Results for the nozzle region of the facility are not 
included, except for a few comments, because of space 
constraints. Including the nozzle region of the facility 
is simply done by including the area changes associated 
with the nozzle. Corrections to the nozzle geometry 
must be included to account for the effect of the large 
boundary layers in the nozzle. 

Expansion Tube Simulation 

A unique feature of the expansion tube among dif- 
ferent types of pulse facilities is that it theoretically 
avoids stagnating the flow within the facility and, 
therefore, avoids any excitation and dissociation of the 
test gas. If this can be done, experimental conditions 


matching the freestream conditions seen by a flight 
vehicle can be produced. 15 Figure 13 depicts the com- 
ponents and ideal operating sequence (using an x-t 
diagram) of the HYPULSE expansion tube which was 
previously the NASA Langley expansion tube . 16 - 17 It 
is composed of the following three sections: a driver, 
an intermediate section (containing the test gas), and 
an acceleration section. For the particular case consid- 
ered here, the driver was 2.44 m (8.01 feet) in length, 
the intermediate section was 7.49 m (24.57 feet), and 
the acceleration tube was 14.62 m (47.96 feet). The 
diameter of the driver is 16.51 cm (6.5 inches) and the 
diameter of the intermediate and acceleration sections 
is 15.24 cm (6 inches). The facility is driven by com- 
pressed helium, and the gas mixtures and pressures 
in the intermediate and acceleration tube are varied 
to achieve the prescribed test conditions. Typical test 
times are from 200 to 600 microseconds. 

As in the previous facilities, the operation of the ex- 
pansion tube begins by rupturing a diaphragm (in this 
facility a double diaphragm) separating the driver gas 
from the test gas. An incident shock wave travels into 
the test gas, compressing it. At the end of the inter- 
mediate tube, a secondary diaphragm is encountered 
and is ruptured. This creates a second incident shock 
and a second expansion both of which travel down the 
acceleration tube. The secondary expansion is prop- 
agating upstream into the test gas but is convected 
downstream in the supersonic flow. The unsteady ex- 
pansion of the test gas creates the high velocity test 
conditions. Data is taken at the end of the acceler- 
ation tube in the gas between the secondary contact 
discontinuity and the tail of the secondary expansion. 
In this ideal operating sequence, no stagnation regions 
are created and, thus, there are no high-temperatures 
to dissociate the test gas. 

Several numerical simulations are presented for the 
HYPULSE facility at a single operating condition. 
One simulation is for the ideal operating conditions 
as explained above and the others investigate the ef- 
fects of a delayed opening of the secondary diaphragm. 
This latter situation is studied because there is clear 
evidence from the pressure transducers positioned near 
the secondary diaphragm that the diaphragm does not 
open instantaneously. Several works have studied the 
effect of the finite opening time on test flow properties 
such as wall pressure, pitot pressure, and velocity. 16,18 
However, the effect of this behavior on the flow qual- 
ity with regard to species concentrations has not been 
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studied in detail. A delay in the diaphragm opening 
causes the primary shock to reflect off the diaphragm 
and create a small stagnation region of hot gas in the 
intermediate tube (the test gas). This will produce dis- 
sociation of the test gas which could impact the quality 
of the flow. The finite-rate chemistry capability in the 
present code allows the amount of dissociated species 
in the test flow created from such non-ideal behavior 
to be estimated. The finite opening time of the di- 
aphragm is simulated by noting the time of arrival of 
the primary shock at the secondary diaphragm and 
then instantaneously removing the diaphragm when a 
specified time has elapsed (i.e. 10 microseconds). No 
viscous or heat-conduction source terms are included 
because the validity of these terms is questionable for 
the high speed, low density flow in the acceleration 
tube. Some indication of the character and size of the 
viscous phenomena which occur inside the accelera- 
tion tube is shown in the axisymmetric computations 
of Jacobs 19 who simulated the HYPULSE facility us- 
ing a perfect gas code. 

For the particular HYPULSE conditions of interest 
here, the driver was filled with unheated helium and 
both the intermediate and acceleration section con- 
tained air. The reported initial driver pressure was 
37.9 MPa (5500 psi), the intermediate section pres- 
sure was 3426 Pa (.497 psi), and the acceleration sec- 
tion pressure was 7.2 Pa (.001 psi). These same val- 
ues were used in the simulations. The initial temper- 
atures used in the simulations were 380K, 292K, and 
292K for the driver, intermediate, and acceleration sec- 
tions, respectively. The temperature of the driver was 
not measured in the experiments but is known to be 
somewhat higher than room temperature because of 
Joule-Thomson and compression heating. 16 A driver 
temperature of 380K was adopted because this tem- 
perature produced primary shock speeds in the inter- 
mediate tube that were close to those measured in the 
experiment. 

Figure 14 contains an x-t diagram of the logarithm 
of density for the flow within the expansion tube with 
a 10 microsecond opening time for the secondary di- 
aphragm. This flow can be compared to the ideal x-t 
diagram in Figure 13. The enlargement of the sec- 
ondary diaphragm region depicts the incident shock 
reflected by the rigid diaphragm and shows that the 
reflected shock is very quickly weakened by the sec- 
ondary rarefaction after the diaphragm is removed. 
However, as seen in the full scale diagram, the re- 


flected shock is still able to cause a significant dis- 
turbance when it reaches the primary interface. Part 
of the wave is transmitted through the interface and 
part is reflected in the form of the shock. This kind 
of interaction is observed and computed in Miller and 
Shinn. 18 

Figure 15 contains a experimental wall pressure 
trace for the conditions mentioned above. 20 Figure 16 
contains computed pressure traces at the test section 
for the ideal case with no delay in the opening of the 
secondary diaphragm, a case with a 10 microsecond 
delay, and a case with a 25 microsecond delay. The 
computed pressure values are slightly higher than the 
experimental values which is expected since no viscous 
correction (no shock attenuation) is included in the 
simulation. All of the computations give similar pres- 
sure traces except for the presence of a shock that has 
reflected off the primary interface for the 25 microsec- 
ond case. For longer opening times (not shown), this 
reflected shock overtakes the test flow and effects the 
test time. These results match the experimental ob- 
servations of Miller and Shinn 18 which were made by 
increasing the secondary diaphragm thickness to ex- 
perimentally increase the opening time. In their work, 
the wall pressures were not greatly effected by the di- 
aphragm thickness, but others quantities such as pitot 
pressure did change. 

A further validation of the computations is made 
by comparing experimental 21 and simulated pressure 
traces from a location 3 inches downstream of the sec- 
ondary diaphragm found in Figures 17 and 18, respec- 
tively. The pressure spike in Figure 17 which occurs 
just after the passage of the secondary shock is caused 
by the finite diaphragm opening. Note that computed 
trace with an instantaneous diaphragm opening has no 
pressure spike. Similar experimental plots for different 
gases and diaphragm thicknesses are found in Miller. 16 
As the diaphragm opening time is increased, the size 
of the disturbance increases. Even with the very sim- 
ple model for the opening time used in this work, the 
qualitative features of the disturbance compare well 
with the experiments. A quantitative comparison of 
the magnitude and the length of time of the distur- 
bance can be used to estimate the size of an “equiv- 
alent opening time” using the current model. In this 
way, the typical opening time present in the experi- 
ments appears to be between 10 and 25 microseconds. 

Figure 19 contains traces of O mass fractions from 
the simulations. The large amount of monatomic oxy- 
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gen immediately behind the secondary shock is caused 
by the very high acceleration gas temperature behind 
the secondary shock. The region of the test gas is 
behind the secondary interface and is labeled in the 
figure. For the ideal case there is very little dissocia- 
tion, as expected. However, for both the 10 and 25 mi- 
crosecond opening times there is a significant amount 
of dissociation in the test gas. The computed traces 
of NO mass fractions, which are not shown here, show 
that this molecule is present in significant amounts, as 
well. 

The high mass fractions of dissociated species in the 
test flow originate at the secondary diaphragm from 
dissociation behind the reflected primary shock. The 
high temperatures and pressures behind the reflected 
shock bring the gas to the equilibrium conditions al- 
most immediately. These high mass fractions persist 
until the test region because the secondary expansion 
freezes the chemical composition before full recombi- 
nation can occur. The simulations indicate that be- 
cause the expanded test flow originates from gas im- 
mediately adjacent to the secondary diaphragm, the 
area affected by a reflected shock does not have to be 
large to significantly impact the test conditions pro- 
vided by the facility. 

The computations presented here were calculated 
with a mesh of 1400 points. Nearly identical results 
were obtained with 700 point and 2800 point grids. 
Solution times for the 1400 point mesh were approxi- 
mately 20 minutes on the NASA Ames Cray Y-MP. 

It is recognized that much more could be said about 
the expansion tube simulations and better modeling 
could be made in several areas: a more sophisticated 
diaphragm opening model could be used (the rupture 
process is surely a multi-dimensional phenomenon), 
boundary layer effects could be taken into account, 
and finite-rate vibrational relaxation could be incor- 
porated. However, the current results agree with ex- 
periments closely enough to suggest that significant 
concentrations of dissociated species may be created 
in the test gas for even relatively small finite opening 
times of the secondary diaphragm and suggest that 
further investigations should be carried out. 

Conclusion 

A quasi-one-dimensionai methodology for numeri- 
cally simulating the flow inside high enthalpy pulse 
facilities was presented. The numerical formulation 
solves the Euler equations with 14 species equations 


with fully-coupled finite-rate chemistry. A moving 
mesh and tracking of gas interfaces were included to 
overcome some of the numerical difficulties associated 
with these types of flows. The value of simulating the 
flow inside a complete facility using the current nu- 
merical approach was demonstrated by computing the 
flows through three different types of pulse facilities. A 
simulation of a helium driven shock tube showed that 
computations can be used to predict the off-tailored 
behavior of shock tubes and tunnels. Computations 
of the flow through the NASA Ames 16-Inch com- 
bustion driven shock tunnel showed the influence of 
non-uniformities in the driver section on the reservoir 
conditions. Finally, the effect of finite secondary di- 
aphragm opening times on the chemical composition 
of the test flow in the HYPULSE expansion tube were 
investigated. 
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Fig, 1 Shock tube temperature profile; initial pressure Fig. 2 Shock tube temperature profile; initial pressure 
ratio= 10, equally-spaced grid. ratio= 10,000, equally-spaced grid. 
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Fig. 3 Shock tube temperature profile; inital pressure 
ratio= 10,000, initial grid clustering. 



Fig. 6a Pressure trace comparison between experi- 
ment and computation for 300 torr case. 



Fig. 4 Shock tube temperature profile; inital pressure Fig. 6b Pressure trace comparison between experi- 

ratio= 10,000, tracking and initial grid clustering. ment and computation for 100 torr case. 


Flux Computed with TVD Scheme 



Cell Side Translated at Local Fluid Velocity 



Fig. 5 Schematic of finite volume cells near a gas Fig. 6c Pressure trace comparison between experi- 
interface to illustrate the inviscid flux calculations. ment and computation for 20 torr case. 
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Fig. 9 Experimental pressure traces for a closed driver 
shot. 



Fig. 12 Comparison of experimental and computed 
pressure traces at DR3 (near diaphragm). 
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Figure 13: HYPULSE Expansion Tube. 
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Fig. 17 Experimental pressure trace 3.0 inches down- 
stream of the secondary diaphragm. 
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VALIDATION OF THREE-TEMPERATURE NOZZLE FLOW CODE NOZ3T 


Seung-Ho Lee, Eloret Institute, Palo Alto, CA 

and 

Chul Park, NASA Ames Research Center, Moffett Field, CA 

Abstract 

A computer code N0Z3T which calculates one-dimensional flows of partially dissoci- 
ated and ionized gases in an expanding nozzle developed earlier is validated in this work 
by comparing with three sets of existing experimental data. The code accounts for the 
differences among the translational-rotational, vibrational, and electron-electronic temper- 
atures. Certain reaction rates are assumed to be controlled by vibrational and electron- 
electronic temperatures. The calculations are compared with: the electron temperature 
and number density data taken in a shock tunnel at CALSPAN during the 1960s, the 
mass-spectrometric data taken in an arc-jet wind tunnel at AEDC during the 19G0s, and 
the spectroscopic data taken in an arc-jet wind tunnel at Ames Research Center recently. 
The results show that the three sets of data can be numerically Reproduced by using the 
code. Typically, vibrational temperature is slightly higher than translational-rotational 
temperature, and electron-electronic temperature is considerably higher than vibrational 
temperature. Atomic oxygen concentration at the end of expansion is significantly smaller 
than that calculated by the conventional one-temperature model. 

Method 

A computer code NOZ3T has been developed in Ref. 1. This code enables one- 
dimensional calculation of a flow expanding through a convergent-divergent nozzle in the 
dissociated and ionized regime. It accounts for three temperatures, that is, translational- 
rotational, vibrational, and electron-electronic. The interactions among these three tem- 
peratures are calculated generally using the model described in Ref. 2. However, three 
different innovations have been made. First, the forward rates for the reactions 

N 2 + 0 -► NO + N 
NO + O — ► 0 2 + N 

are assumed to be determined by the electron-electronic temperatures. This is because 
these reaction rates are controlled mostly by the concentration of electronically-excited 
O-atoms. Second, radiative loss is accounted for. Third, reactions involving N0 2 and N 2 0 
are accounted for. 


Results 

Ames Arc-Jet Data - An experiment has been conducted recently at NASA Ames 
Research Center to determine temperature of the freestream flow in the test section of an 
arc-jet wind tunnel. The tunnel was operated at the stagnation chamber pressure of 2.5 
atm, and enthalpy of about 45 MJ/kg. Emissions from the Gamma and Delta bands of 
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NO molecules are measured spectrographically in the test section. The results showed that 
the rotational temperature is about 600 K, the vibrational temperature is about 900 K, 
and the electronic temperature is about 7000 K for the experimental condition. The code 
NOZ3T is run to numerically reproduce these measured values. By dividing the vibrational 
relaxation time value of Millikan and White by 10 7 , the observed values are reproduced 
numerically by the NOZ3T code. 

CALSPAN shock tunnel data - The temperature and number density of electrons had 
been measured along an expanding nozzle of a shock tunnel at CALSPAN (Refs. 3 to 5). 
The experimental data are reproduced by NOZ3T. The comparison between the measured 
and the calcualted values is given in Figs. 1 to 3. As seen in these figures, NOZ3T can 
reproduce the experimental data fairly well. Here again, vibrational temperature is higher 
than the translational temperature, and electron-electronic temperature is higher than 
vibrational temperature. 

AEDC arc-jet data - Concentration of 0, N, N 2 , 0 2 , and NO had been measured in 
the test section of an arc-jet wind tunnel at AEDC (Refs. 6 and 7). Using the NOZ3T 
code, the measured data have been numerically reproduced. Comparison between the 
experimental and the theoretical results is given in Figs. 3 to 6. In these figures, the solid 
curves represent the theoretical values obtained using the conventional one-temperature 
model. As seen in these figures, one-temperature model is not able to reproduce the 
measured data. The present three-temperature model can reproduce the data fairly well. 
It is important to note that measured concentrations of atomic oxygen is less than those 
calculated by the one-temperature model. The three-temperature model correctly predict 
this trend. 


Conclusions 

Using NOZ3T code, it is possible to numerically reproduce the three existing ex- 
perimental data in which the nonequilibrium thermochemical state of the gas has been 
determined at the exit of a hypersonic nozzle. Generally, vibrational temperature is higher 
than translational-rotational temperature, and electron-electronic temperature is higher 
than vibrational temperature. Concentration of atomic oxygen is less than that calculated 
by the conventional one-temperature model in typical environments. 
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FIGURE CAPTIONS 

Figure 1. Comparison between the experimental data and the present calculation for shock 
tunnel flow with nitrogen; settling chamber conditions: p = 17.1 atm, T = 7200 K (Ref. 

4 ) 

Figure 2. Comparison between the experimental data and the present calculation for shock 
tunnel flow with oxygen; settling chamber conditions: p = 25 atm, T = 4950 Iv (Ref. 5) 

Figure 3. Comparison between the experimental data and the present calculation for shock 
tunnel flow with air; settling chamber conditions: p = 25 atm, T = 7200 Iv (Ref. 3) 

Figure 4. Comparison between the experimental data, and the present calculation for 
arc-jet flow with air; nitrogen and oxygen atoms (Ref. 6 and 7). 

Figure 5. Comparison between the experimental data and the present calculation for 
arc-jet flow with air; nitric oxide (Ref. 6 and 7). 

Figure 6. Comparison between the experimental data and the present calculation for 
arc-jet. flow with air; nitrogen and oxygen molecules (Ref. 6 and 7). 
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Abstract 

The flow field around a proposed lunar return aerobrake is examined computationally assuming viscous, 
laminar flow and utilizing an effective y approach to incorporate real-gas effects. The flow fields for diree cases 
are calculated with both an axisymmetric and a three-dimensional formulation. The three vehicle configurations 
have braking panels extended 50°, 55°, and 60° to the inflow (zero lift configurations). The resulting 
axisymmetric and 3-D flow fields are examined, and it is shown that despite complexities in the 3-D flow fie d, 
the total drag coefficients calculated from modified axisymmetric results are very close to those obtained from the 
3-D solutions. The aerobrake is shown to achieve total drag coefficients as high as 8.4 for a vehicle with panels 

deflected to 60°. 


Introduction 

An aerobrake uses aerodynamic forces, in 
place of rocket propulsion, to decelerate and 
maneuver, thereby achieving orbital modifications 
or entry into the atmosphere. The aerobrake 
concept examined here is intended for lunar return 
missions (11 km/sec entry velocity). It is designed 
so that, in a stowed configuration, it can be 
launched into space using the proposed National 
Launch System (NLS) or the Space Shuttle which 
has a payload bay capacity of 4.57m x 18.28m. 

Nomenclature 

A reference area, rcd 2 /4, m 2 

Cp total drag coefficient, D/q„-A 

D drag.N 

d diameter, m 

L lift, N 

M mass, kg 

q*, dynamic pressure, N/m 2 

r radius, m 

5 panel deflection angle 

y specific heat ratio, Cp/C v 


* Research Scientist, Member AIAA. Mailing Address: 
NASA Ames Research Center, MS 230-2, MoffeU 
Field, CA 94035 

t Research Scientist, Member AIAA. 

* Research Scientist, Associate Fellow, AIAA. 

This paper is declared a work of the U.S. Government 
and is not subject to copyright protection in .the 
United States. 


Vehicle Design 

An aerobrake is required to generate large 
amounts of drag. Therefore, this initial design 
study concentrates on the prediction of the drag 
capabilities of the proposed vehicle. The aerobrake 
vehicle, shown in Fig. 1, is composed of an hemi- 
ellipsoidal nose (forebody), a nearly cylindrical 
main-body, and four large, flat panels. The panels 
are used mainly as braking surfaces; however they 
are independently adjustable so that lift can be 
generated for maneuvering within the atmosphere to 
alleviate deceleration loads and heating. Desirable 
L/Ds for lunar return are 0.25 to 0.3. The panels 
can be stowed in the launch vehicle with zero 
deflection, and the maximum width and length of 
the panels is limited by the size of the launch 
vehicle's payload bay. Initially, a vehicle with four 
braking panels is examined. The optimum number 
of panels has yet to be determined, and future work 
may address this design feature. The panels can be 
deployed at variable angles to the freestream 
depending on the amount of drag required. Three 
panel deflection angles are examined in order to 
study how the panel angle affects the total drag 
coefficient 

The radius of the main-body increases from 
the forebody to the end of the main body. 
However, parts of the main body are flattened to 
provide a long hinge line for attachment of the 
panels. Tapering of the lower portion of the panels 
and the flattening of the main body allow for a long 
hinge line so that loads on the panels can be 
distributed along the hinge line without unduly 
reducing the volume of the vehicle. 
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The panels are made large to yield a low 
ballistic coefficient, M/Cq-A. This permits the 
vehicle to decelerate at higher altitudes thereby 
reducing deceleration, structural loads and heating 
rates. The hemi-ellipsoidal nose shape was chosen 
to maximize drag and minimize heating to the 
forebody region. The payload capacity of the 
aerobrake is estimated to be 10,000 kg, and it has a 
total volume of about 150 m 3 . The empty weight 
of the aerobrake is assumed to be 6500 kg, and the 
entry weight is 16,500 kg. 

Approach 

For initial design purposes, a general 
understanding of the over all flow field 
characteristics and the drag producing capabilities of 
the aerobrake is required. For this initial study, 
all four panels are maintained at the same deflection 
angle (zero lift configuration), and three panel 
orientations are examined: 5=50°, 55°, and 60°, 
The range of panel angles was chosen based on past 
drag brake studies 2 which showed that for similar 
vehicles with no gaps in the braking surfaces, these 
deflection angles yielded the highest drag while 
maintaining good flow quality and relatively low 
heating rates. 

The freestream (flight) conditions used for 
this study (see Table 1) correspond to the peak 
pressure point along an aerobrake trajectory. A 
ballistic coefficient, f}=250, and a maximum lift to 
drag ration, L/D=0.25, were assumed to calculate a 
lifting trajectory for the aerobrake. For the fi and 
entry mass assumed, the corresponding Cd is 
approximately 4. The L/D for this trajectory is not 
held constant, and the maximum g-loading for this 
trajectory is 8.3gs. 

Note that for all calculations requiring a 
reference area, the area is based on the largest main 
body diameter (d=4.572 m). 

The flow fields for these cases are 
calculated assuming a laminar boundary layer and an 
ideal gas having an effective constant specific heat 
ratio. It is also assumed that the pressure in the 
base flow region is very low and will not affect that 
calculation of total drag coefficient The panels 
create compression comers that, at the hypersonic 
speeds examined, will cause the flow to massively 
separate. Because of the separation, Newtonian 
flow or other inviscid methods are not applicable to 
this problem, and full Navier-Stokes solvers are 
required. 


A large portion of the aerobrake trajectory 
is dominated by real-gas effects. These effects are 
approximated by incorporating an effective value of 
the specific heat ratio, y, in the code. This method 
has been shown to be suitable for the calculation of 
drag and lift coefficients of blunt bodies in real-gas 
flows.* » 2 At the chosen flight conditions the 

stagnation pressure is nearly 0.4 atm; therefore the 
flow is assumed to be in chemical and thermal 
equilibrium. An effective specific heat ratio, y, of 
1.21 was used 3 based on the flight conditions. 

Because the panels are flat and portions of 
the main-body are flattened, the aerobrake geometry 
is three dimensional. However, it is very nearly 
axisymmetric, and therefore, an axisymmetric 
geometry is used as a first approximation to reduce 
the time required for the flow field computations. 
The axisymmetric geometry is created by generating 
a body of revolution using a panel symmetry line. 
The surface pressure field for this axisymmetric 
geometry is then computed. The calculation of drag 
is a function of the vehicle geometry as well as the 
surface pressure field. Therefore, to better estimate 
the drag on the 3-D vehicle using the axisymmetric 
results, the axisymmetric pressure field is applied to 
the actual 3-D geometry. In this way, the flatness 
of the panels, the flattened portions of the main 
body, and the gaps between the panels are taken 
into account The drag coefficient that is calculated 
by integrating the axisymmetric pressure field over 
the 3-D geometry is then referred to as a "modified" 
axisymmetric value. 

In addition, solutions are also computed 
using the actual three-dimensional geometry. This 
geometry contains one octant of symmetry (45° 
from the middle of a panel to the middle of a gap 
region); therefore, only one eighth of the full 
configuration is needed for the three-dimensional 
calculations. 

Numerical Methods 

Two different CFD codes were used to 
calculate the axisymmetric flow around the 
aerobrake. The first code is an axisymmetric solver 
developed by Candler .4 This is a finite-volume 
code that uses Steger- Warming flux-vector splitting 
and Gauss-Seidel line relaxation to solve the 
Navier-Stokes equations. This code has been used 
previously for similar studies of hypersonic drag 
brakes. 2 

The second code that was used is a 3-D 
solver, FL3D, developed by Venkatapathy. 5 This is 
a finite difference code that uses a three- 
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dimensional, implicit time marching Navier- 
Stokes solution procedure. This code has 
previously been applied to axisymmetric and three- 
dimensional nozzle-plume and blunt body flows 5 . 
The numerical method is a LDU-ADI scheme with 
Roe's averaging and MUSCL differencing. High 
order spatial accuracy is achieved by constructing 
MUSCL extrapolation of flow variables with a 
differentiable limiter. Roe averaged upwind fluxes 
are constructed and the flux -difference from the 
higher order fluxes are used in this formulation. 
The LDU-ADI algorithm is a diagonal algorithm 
requiring minimal CPU per iteration and is 
applicable to steady flows. Both the codes are ideal- 
gas solvers which incorporate an effective value of 

Y- 

Since no experimental data exist for the 
present flow Held, solutions from the two codes 
were compared against one another to assure that 
the flow fields are calculated consistently. The flow 
field solutions for an axisymmetric aerobrake with 
panels at 50° were used for the comparison. The 
general flow structures, including the size of the 
separation bubble caused by the panels, compared 
well, and the differences in calculated surface 
pressures were less than 1%. The size of the 
separation region is important because it affects the 
efficiency of the braking panels. The pressure 
within a separated region remains essentially 
constant even when part of that region is a 
compression surface such as the base of a panel. 
Therefore, the portion of the panel contained within 
the separated region contributes little to the total 
drag. 

For each set of panel angles, an 
axisymmetric and a 3-D solution were calculated. 
Although the 3-D solver could be used to calculate 
the axisymmetric solutions as well as the 3-D 
solutions, the axisymmetric solver was used 
because of its speed and convenience. Each 
axisymmetric solution was started from the same 
initial conditions and required approximately 10-15 
minutes of CPU time on a Cray Y-MP. The first 
3-D solution, the 50° panel case, was started from a 
blunt cone solution and required approximately 1.75 
hours of CPU time. Subsequent 3-D solutions 
were started from the most similar preexisting 3-D 
solution. For example, the solution from the 50° 
panel case was used as the initial conditions for the 
55° panel case. This reduced the CPU time to 
approximately 45 minutes. 

The grids used are 63x47x61 for the 3-D 
calculations and 63x61 for the axisymmetric 
calculations. A representative axisymmetric grid is 
shown in Fig. 2. For the axisymmetric 


calculations, a solution-adapted 6 grid technique was 
used to refine the major flow features. Although 3- 
D grid adaptation is possible with the same 
procedure, it was concluded that solution adaptation 
of the 3-D grid was not necessary. This conclusion 
is based on the axisymmetric results which showed 
that although the flow features became more 
distinct, the calculated drag coefficient changed very 
little from the non-adapted cases. For the 3-D 
calculations, symmetry is taken advantage of and 
only one octant (45°) of the full 3-D geometry is 
computed. The surface grid for the 3-D 
calculations, shown in Fig. 3, consists of half of a 
panel and half of a gap region. The forebody of the 
aerobrake is truly axisymmetric at the zero angle of 
attack cases considered here and is not influenced by 
the downstream flow. Therefore there is no need to 
compute the nose region using a 3-D formulation. 
The portion of the flow field that is actually 
computed by the 3-D solver begins at the end of the 
forebody, and the inflow is fixed to the values from 
the axisymmetric solution. This reduces the grid 
size for the 3-D calculations to 46x47x61. 

The outflow boundary conditions used for 
the gap region of the 3-D geometry are of concern 
because a portion of this region may be 
recirculating, subsonic flow. The supersonic 
portions of the outflow are extrapolated from 
interior points. The pressure is fixed across the 
subsonic region to the pressure that exits just as the 
flow becomes supersonic. The remaining flow 
variables are then extrapolated. Fixing the pressure 
across the subsonic region is a fairly good 
assumption because the pressure within separation 
bubbles remains relatively constant. This is indeed 
the case for the portion of the separation region 
which affects the aerobrake surface pressure. 

Results and Analysis 

Colored pressure contours for a 
representative axisymmetric solution at 8=50° are 
shown in Fig. 4. The pressure contours reveal that 
the peak pressure on the panels occurs in the region 
where the forebody shock interacts with the shock 
formed off the panel. The separation and 
reattachment regions are also indicated. Notice that 
the separation region is quite large. Consequently 
the pressure on the panel surface remains low 
within the separation zone and only begins to 
substantially increase downstream of the 
reattachment point 

A comparison of the axisymmetric and 3- 
D calculations of surface pressure is made for an 
aerobrake with panels at 50° (see Fig. 5). 
Although the axisymmetric solution shows a 
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somewhat higher peak pressure on the panel, the 
surface pressure profiles look very similar 
indicating that, at least along a panel centerline of 
the 3-D geometry, the axisymmetric solution 
comes very close to the 3-D case. However, for 
either case it is noteworthy that the multiple shock 
interaction (see Fig. 4) increases the peak pressures 
on the panels by about a factor of 2 over the 
Newtonian impact theory value of nearly 23 kPa. 

The limitations of the axisymmetric 
solution become apparent when numerically 
simulated surface oil flow for the 3-D solution is 
examined (see Fig. 6). The surface pressure on the 
body is indicated by colors; dark blue is the lowest 
pressure, essentially freestream, and white the 
highest The oil flow is shown in yellow. The 
separation on the main body indicated by the 
convergence of oil flow lines is largest in the center 
of the panel regions and decreases only slightly in 
the centers of the gap regions. The patterns of the 
oil flow indicate that a complex, multiple vortex 
structure exists within the separated region. 

Particle traces constrained to a panel center 
plane are shown in yellow in Fig. 7. The particle 
traces indicate the size of the separated region and 
show more clearly than the surface oil flow the 
extent of the separation on the panels. Almost half 
of the panel length is contained within the 
separation region. The effect of the separation on 
the panels is demonstrated by the painted surface 
pressure contours which show a very low pressure 
within the entire separated region. A cut away of 
the outer shock region is shown in this figure as a 
solid surface. The effect of the separation on the 
outer shock surface is indicated by the slight 
bulging of the surface in the region of the 
separation. Because of the panel-gap geometry, the 
shock structure in the flare region is truly three 
dimensional. It retains a conical shape in the gap 
regions and bulges out in the panel regions where it 
reflects off the panels. The shock stand off distance 
on the panels is greatest at the center of the panels. 
As the flow expands around the edge of the panels, 
the shock curves and the stand-off distance 
decreases. This causes the peak pressure on the 
panels to shift off the center line to a position 
approximately midway between the panel center and 
the edge. 

Particle traces constrained to a gap center 
plane are shown in Fig. 8. The particle traces 
terminate at the end of the computational domain. 
The separation in the gap region remains relatively 
large and, as is indicated by the traces, the 
recirculation region extends into the base area where 
the flow field has not been calculated. The flow at 


the outflow boundary of the computational domain 
remains mostly supersonic even within the 
separation, therefore, the extrapolated boundary 
conditions applied at the outflow region are 
adequate. However, the effect of the base flow 
region on the main body is not captured by this 
representation. 

It is assumed that the base flow region has 
very little affect on the total drag that is calculated 
for the vehicle, since the pressure in that region is 
relatively low. However, because of the gaps, the 
base flow may have some effect on the main body, 
and it is not clear what these effects may be. To 
capture the base flow region properly, the flow 
behind the aerobrake needs to be calculated for a 
distance equivalent to approximately four times the 
panel’s vertical height. This corresponds to the 
approximate distance required for the flow to 
reattach (if the reattachment point is not contained 
within the computational domain, the error is not 
eliminated). Increasing the grid size to incorporate 
the necessary base flow region would increase the 
solution time substantially. The major drag 
producing areas of the vehicle are on the panels 
where the shock-shock interaction occurs. It is 
unlikely that including the base flow region in the 
calculation would affect this part of the flowfield. 
Therefore, the increase in solution time needed to 
calculated the base flow region is not warranted in 
this initial study. 

A comparison of the total Cp calculated 
for the modified axisymmetric and 3-D solutions for 
50°, 55°, and 60° flare angles is shown in Fig. 9. 
Here the term "modified" indicates that the 
axisymmetric Cq is computed by first distributing 
the axisymmetric surface pressure field across the 3- 
D geometry, thus taking into account that there are 
gaps and that the panels are flat This configuration 
is then used to calculate the "modified" 
axisymmetric Cp. The most notable feature of this 
comparison is that, although the axisymmetric 
solution is incapable of computing the complex 
flow field revealed by the 3-D solutions, the total 
drag coefficient calculated for the modified- 
axisymmetric and 3-D solutions are very close. 
This suggests that, for this particular panel-gap 
geometry, it is not necessary to compute a full 3-D 
solution in order to calculate drag coefficient. 
However, it is not clear from this study whether the 
relationship between the axisymmetric and 3-D 
results would be maintained if the gap size was 
changed, if the total number of panels was increased 
or decreased, or if the panel angles were different 
from those examined here. For the cases 
considered, the general structure of the flow fields 
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did not change substantially with the changes in 
panel angle. 

The main purpose of the panels is to 
increase the total drag of the vehicle. Figure 10 
shows how the (modified) drag coefficient changes 
along the vehicle's length for all three panel 
deflection angles. The base of the panels is at 
x=11.23 id. Because the vehicle’s main body is 
cylindrical for these axisymmetric cases, the main 
body does not contribute anything to the total drag, 
and the Cd curve is flat in this region. The curve 
does not begin to rise until x=15m. This indicates 
that approximately (15m-11.23m=) 3.77m or 
53.5% of the panels are contained within the 
separation region and contribute little to the total 
drag. 

The total drag coefficient of the vehicle 
without the panels is Cd=0 886. The vehicle with 
panels at deflection angles of 50°, 55°, 60° achieves 
total drag coefficients of 5.84, 7.24, and 8.39, 
respectively. This corresponds to an increase in 
total Cq between 560% and 847% above the no- 
panel configuration. Also note that the total Cd 
increases 24% from the 50° deflection case to the 
55° deflection case, and the increase from the 55° to 
60° case is 16%. The resulting changes in the 
vehicle's ballistic coefficient are shown in Fig. 11. 

The amount of the total Cd increase can 
be tailored by changing the panels' lengths and 
deflection angles. The panels’ length would be fixed 
for a specific mission category (i.e. lunar return or 
mars return) and the total Cd variation needed 
during a specific mission would be achieved by 
changing the deflection angles. 

The effects of turbulence were not 
examined in this study; however a turbulent 
boundary layer tends to decrease the separation 
bubble size thus causing the total drag to increase. 
Therefore, the laminar calculations of drag 
coefficients are conservative estimates. 

Conclusions 

The 3-D calculations of a proposed lunar 
return aerobrake show that a major portion of the 
main body as well as the braking surfaces are 
contained within a large, complex separated flow 
region. Despite the complexity of the 3-D flow 
field, modified axisymmetric solutions predicted 
total drag coefficients very close to those that were 
calculated from the 3-D solutions for aerobrakes 
with panels at 50°, 55° and 60° angles. 


It was shown that, for the panel deflection 
angles and panel length examined here, the total 
drag coefficient could be increased to as much as 8.4 
times the Cd of the same vehicle without panels. 
The increase in total Cd that is required for a given 
mission can be achieved by changing the lengths 
and deflection angles of the panels. It was also 
shown that Newtonian impact theory predicted a 
peak panel pressure one half of that calculated by 
the present methods showing that the Navier-Stokes 
method w as indee d necessary. 

Questions were raised concerning the 
possible effects of the base flow region and also 
whether the agreement between the axisymmetric 
and 3-D calculations would be maintained for 
similar aerobrakes with different panel 
configurations. These questions will be examined 
in future work 
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altitude 

56.6 km 

velocity 

9082 m/s 

density 

4.751x 10~4 kg/m3 

pressure 

35.0 N/m2 

temperature 

256.6 K 

wall temp. 

2700 K 


Table 1 Flight altitude, velocity, and ambient 
conditions 








Figure 8 Calculated surface oil flow and surface pressure for 5=50‘ 
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Figure 9 Total drag coefficient comparisons Figure 10 Variation of drag coefficient down the 

between axisymmetric and 3-D calculation length of the vehicle for three panel deflection angles 



Figure 1 1 Change in ballistic coefficient due to panel deflection 
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Abstract 


An initial numerical/experimental investigation has 
been done to better understand multi-dimensional flow 
phenomena inside pulse facilities. Time-dependent 
quasi-one-dimensional and axisymmetric numerical 
simulations of complete shock tube flow are com- 
pared with experimental pressure traces recorded at 
the NASA Ames electric-arc driven shock tube facil- 
ity (from cold driver shots). Of particular interest is 
the interaction between the reflected shock wave and 
the boundary layer. Evidence of the shock bifurcation 
caused by this interaction is clearly seen in the present 
experimental data. The axisymmetric simulations re- 
produce this phenomenon and demonstrate how this 
interaction can provide a mechanism for driver gas to 
contaminate the stagnation region. The simulations 
incorporate finite-rate chemistry, a moving mesh and 
laminar viscosity. 

I. Introduction 

The duration of the test time provided by a shock 
tube or shock tunnel is always less than the theoret- 
ical value. In fact, it is not uncommon for the test 
time to be only a fraction of what might be ideally 
expected. The usual reason for this decrease is the 
early arrival (or contamination) of driver gas in the 
stagnation region, although there may also be other 
causes such as the arrival of the rarefaction wave re- 
flected off the driver end of the shock tube. Since 
accurately knowing the experimental conditions pro- 
vided by such facilities is important, many investiga- 
tions have been carried out to understand and quan- 
tify the physical mechanisms which cause shortened 
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test times. Several of these mechanisms are schemat- 
ically presented in Figure 1 and include deformation 
of the contact discontinuity caused by the diaphragm 
rupture process, 1 mass transfer of the driven gas into 
the boundary layer, 2 contact discontinuity diffusion 
and instabilities, 3-5 and shock/boundary layer inter- 
action after reflection of the incident shock off the end 
wall, 6-18 Non-ideal performance may also be caused 
by the shock tube/tunnel geometry. For example, 
Jacobs 19 and Cambier et a/. 20 have numerically pre- 
dicted the formation of a vortex created during the 
reflection of the incident shock off an end wall which 
includes a nozzle entrance. This vortical flow may en- 
hance the mixing of the driver and driven gases in 
shock tunnels. 

It is not clear which of the mechanisms mentioned 
above limits test times most and it may be that the 
dominant mechanism varies with the experimental fa- 
cility or the run conditions. There is a large amount 
of evidence, however, that suggests that the reflected- 
shock/boundary layer interaction is often a major con- 
tributor to the contamination process. It is well estab- 
lished that under many conditions the reflected shock 
will interact with the boundary layer causing it to bi- 
furcate near the wall (see Figure 2). Mark 6 developed 
a model to predict the characteristics of the bifurca- 
tion and the conditions under which it will occur. He 
showed that the flow in the energy deficient bound- 
ary layer has a stagnation pressure that is less than 
the stagnation pressure behind the normal reflected 
shock and is prevented from passing under the re- 
flected shock. Instead, it separates and collects in a 
bubble of gas next to the wall. Davies 8,9 used Mark’s 
model to show that the bifurcated shock structure pro- 
vides a mechanism for driver gas contamination of the 
stagnation region. It can be seen in Figure 2 that the 
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gas which encounters the oblique shocks at the foot 
of the bifurcated shock retains a higher velocity than 
the gas which encounters the normal shock. The re- 
sult is a jet of gas near the wall. Clear experimental 
evidence of this wall jet can be seen in the recent color 
schlieren photographs of Kleine et a/. 12 which include 
features such as the rolling up of the wall jet as it en- 
counters the end wall. Davies showed that when the 
bifurcated shock structure of Figure 2 interacts with 
the contact discontinuity, cold driver gas is propelled 
by this jet toward the end wall and into the driven gas. 
He and others such as Bull and Edwards 10 have done 
experiments which measured the time of arrival of the 
cold driver gas, Davies by heat-transfer gages and Bull 
and Edwards by “tagging” the driver gas with an in- 
frared active gas. Both found that measured values 
of the arrival time of driver gas in the stagnation re- 
gion were in agreement with the analytical theory of 
Davies. An interesting experimental work by Ligong 
et a/. 14 attempted to reduce driver gas contamination 
using boundary layer suction to eliminate the wall jet. 

In the absence of optical data, the presence of shock 
bifurcation can be inferred from side-wall pressure 
measurements. Experiments show that the passage of 
the reflected shock is not accompanied by a sharp rise 
but by a rise of about one half the expected value fol- 
lowed by a small plateau or drop in pressure. The post- 
reflection pressure is then achieved after a pressure 
overshoot (see Figure 3). Sanderson 15 used Mark’s 
model with a small modification to explain how the 
bifurcated foot causes these pressure trace features. 
The presence of reflected-shock/boundary layer inter- 
action can also be deduced by noting the change of the 
reflected shock speed compared to theory. 6, 16,17 

The description of reflected-shock/boundary layer 
interaction presented above is adequate to explain 
much of what is observed experimentally in shock 
tubes. However, more complicated flow structures 
such as a pseudo-shock (or shock train) can develop 
when the shock/boundary layer interaction is strong. 
Matsuo et a/., 16 Strehlow and Cohen, 17 and Brossard 
and Charpentier 18 all show schlieren photographs of 
the formation of multiple shocks after the reflection of 
the incident shock. The effect of these multiple shocks 
on driver gas contamination has not been studied. 

There have been very few multi-dimensional sim- 
ulations of shock tubes and shock tunnels that re- 
solve wall viscous effects and none of them have looked 
at the contamination of the stagnation region with 


driver gas through the wall jet mechanism proposed 
by Davies. The reason for this is surely related to 
the small time steps caused by the fine grid needed 
to resolve boundary layers and the resulting high cost 
of propagating the flow features down the length of a 
shock tube at those small step sizes. Jacobs 21 has sim- 
ulated expansion tube flow of a perfect gas including 
the boundary layer. Cambier et al 20 looked at sev- 
eral isolated multi-dimensional phenomena including 
the diaphragm rupture process. Several other works 
have simulated reflected-shock/boundary layer inter- 
action by considering only the end wall region of a 
shock tube and assuming a boundary layer profile for 
the inflow. 12,22 

The approach taken in this work is to develop a nu- 
merical methodology for time-dependent axisymmetric 
simulations of the flow inside the NASA Ames electric- 
arc driven shock tube (cold driver shots without the 
arc-driver) starting from diaphragm rupture. The nu- 
merical scheme is based on experience gained from 
time-dependent quasi-one-dimensional simulations for 
pulse facilities. 23 The high cost of resolving the bound- 
ary layers and simulating the flow through the entire 
facility is minimized by an efficient use of grid points, 
highly vectorized code, and an implicit treatment of 
the viscous terms. Experimental data is obtained from 
the facility for comparison to the simulations. Starting 
the simulations with the diaphragm rupture allows the 
position of the contact discontinuity and the bound- 
ary layer development to be computed. This, in turn, 
makes it possible for phenomena such as the reflected 
shock/boundary layer interaction and the reflected- 
shock/contact discontinuity interaction to be investi- 
gated numerically. Of particular interest is the ability 
of the simulations to predict the driver gas contami- 
nation by the wall jet mechanism proposed by Davies. 
In this work, the contact discontinuity is assumed ini- 
tially to be planar and the boundary layer is treated 
as laminar. Comparisons of the axisymmetric simu- 
lations with quasi-one-dimensional solutions help il- 
lustrate the pressure trace features captured by the 
axisymmetric solutions. 

II. Experimental Facility 

The NASA Ames electric-arc driven shock tube 
facility. 24,25 has several possible configurations and has 
a large hypervelocity operating range using its arc- 
driver. However, this work only considers experiments 
using a cold helium driver with nitrogen in the driven 
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section. These shots were made with a cylindrical 
driver of .76 m (2.5 feet) in length and 10 cm (3.93 
inches) in diameter (see Figure 4). The driven section 
was 4.02 m (13.17 feet) long with the same diameter 
as the driver. A single self-break diaphragm separated 
the driver and driven gases. 

For the current experiments the instrumentation 
consisted of pressure transducers flush mounted at two 
positions on the shock tube walls. One location was 
.0254 m from the driven tube end wall and the other 
was .445 m from the end wall. The transducers were 
PCB Piezotronics, Inc. Model 113A21 with a circu- 
lar surface area .218 inches in diameter. They have 
a rise time of 1 microsecond. Heat- transfer data has 
also been collected at the same instrument locations 
and will be used in future work. 

m. Numerical Method and Gas Model 

The gas dynamic equations for both the quasi- 
one-dimensional and the axisymmetric simulations are 
solved by using an explicit finite-volume form of 
the Harten-Yee upwind TVD scheme 26 which gives 
second-order spatial and temporal accuracy. The gas 
model includes the three major species present in the 
shock tube for the present experiments (Y 2 , N t and 
He) and accounts for finite-rate chemical processes. 
A separate equation for vibrational energy is included 
so that vibrational nonequilibrium effects can be as- 
sessed. The present work, however, enforces thermal 
equilibrium. The full Navier-Stokes viscous terms are 
included in the axisymmetric simulations. 

In the quasi-one-dimensional computations, a track- 
ing scheme is used for the contact discontinuity to 
avoid the numerical difficulties associated with resolv- 
ing this feature as it moves large distances. 23 The 
tracking scheme is not used in the axisymmetric sim- 
ulations so that mixing of the the driver and driven 
gases can occur during the interaction of the reflected 
shock with the contact discontinuity. Instead, the ax- 
isymmetric simulations cluster grid points at the con- 
tact discontinuity to minimize numerical diffusion and 
convect this clustered grid with the gas interface as it 
travels down the driven tube. This approach has the 
additional benefit of compressing all of the cells asso- 
ciated with the driven tube into the end wall region of 
the shock tube as the driven gas is compressed thereby 
providing a fine axial mesh during the shock reflection. 

Because time accuracy is required, the global time 
step is dictated by the smallest required time step in 


the computational domain. The solutions for both the 
one-dimensional and axisymmetric solutions are ad- 
vanced at a CFL number less than 1 based on the 
inviscid gas dynamics. In the boundary layer of the 
axisymmetric simulations, this CFL number can be 
significantly larger than the stability bound required 
by the viscous terms. To avoid the more limiting time 
step dictated by the viscous terms, they are treated 
implicitly. This requires a block tri-diagonal matrix 
inversion along each line of cells normal to the wall. 
The cost of this inversion is more than offset by the 
greater allowable time step. It is believed that the time 
accuracy of the solution is not significantly effected by 
this approach and without it the simulation becomes 
impractical. The source terms representing the finite- 
rate chemical kinetics and vibrational relaxation are 
also treated implicitly. This implicit treatment reduces 
the formal time accuracy to first order. 

IV. Results 

All of the numerical work in this paper is compared 
with a single experimental condition in the NASA 
Ames electric-arc driven shock tube. The initial con- 
dition for the experiment was 3.96 MPa (575 psi) cold 
helium in the driver and 100 torr nitrogen in the driven 
section. This resulted in an incident shock speed of 
approximately 1490 m/sec. The temperature in the 
driver was not measured and so is not known precisely. 

To illustrate the general character of the flow in the 
NASA Ames electric-arc driven shock tube with the 
cold helium driver, a computed x-t diagram of tem- 
perature contours from a one-dimensional simulation 
is given in Figure 5. The x-t diagram is for a flow at 
a higher shock speed than the conditions mentioned 
above, but the flow is qualitatively similar. This figure 
clearly shows the time- varying positions of the incident 
shock, contact discontinuity, and rarefaction. Included 
in Figure 5 is an enlargement of the end wall region 
of the driven section showing the reflected shock being 
partially transmitted through the contact discontinu- 
ity and partially reflected back toward the end wall. 
This reflection occurs several times and each successive 
interaction further slows the interface down until it is 
brought to rest. This type of off-tailored interaction 
(with reflected shocks) is referred to as an over-driven 
case. Note that a vertical line of data from this plot 
yields a trace over time at a given location. 

Figure 6a shows experimental and computed one- 
dimensional sidewall pressure traces .0254 m upstream 


3 


of the driven section end wall. There are two one- 
dimensional simulations included in the figure to show 
the sensitivity of the solutions to the initial conditions. 
Each of the simulations used 350 points with half in 
the driver and half in the driven section. Both have the 
same the initial conditions except that the driver tem- 
perature is 300 K in one case and 330 K in the other. 
A driver pressure of 3.45 MPa (500 psi) was chosen for 
both simulations to produce approximately the same 
shock speed as measured in the experiment. The re- 
quired driver pressure for the computations is less than 
that used in the experiments because there are losses 
in the experiment due to the diaphragm rupture and 
viscous effects. The shock speeds for the 300 K and 
330 K driver cases are 1495 M/sec and 1529 m/sec, 
respectively. The driver at 300 K produces a smaller 
incident shock speed compared the driver at 330 K 
and nearly matches the measured experimental shock 
speed. The driver at 300 K also slows the arrival of 
the rarefaction off the opposite end of the facility by 
reducing the speed of sound in the driver. Note that 
a relatively small change in the initial conditions has 
a significant effect on the pressure traces. Because 
the exact time of the experimental diaphragm rupture 
is not known, the experiments and computations are 
temporally aligned at the shock arrival time at the lo- 
cation .445 m from the driven tube end wall. 

Several of the major events are labeled in the Fig- 
ure 6a. These include the passage of the incident 
and reflected shocks, the arrival of the rarefaction, 
and the presence of reflected waves off the contact 
discontinuity. The general changes in the pressure 
traces due to off-tailored interactions with the con- 
tact discontinuity are captured fairly well with the 
one-dimensional approach. As expected, the experi- 
mentally measured features associated with the con- 
tact discontinuity (such as waves which have inter- 
acted with it) are smeared because the interface is 
non-planar and diffuse while the same features in the 
simulations are sharp because the interface is tracked. 
The overshoot in the experimental trace just after the 
reflected shock is from the bifurcation of the reflected 
shock as explained by Sanderson. 15 

Figure 6b contains the pressure traces from the ex- 
periment and the two one-dimensional simulations for 
the location .445 m upstream of the driven tube end 
wall. The main feature to consider in this figure is 
the arrival time of the reflected shock at this station. 
Both one-dimensional simulations significantly under 


predict the reflected shock speed compared to the ex- 
periment. The disagreement is expected because the 
reflected shock is significantly influenced by the shock/ 
boundary layer interaction on the shock tube walls 
and because the one-dimensional simulations cannot 
account for this phenomenon. Figure 6b also shows 
that the reflected shock speed is more strongly affected 
by the shock tube initial conditions than the incident 
shock (the 300 K and 330 K cases are quite different). 

Axisymmetric simulations were then performed for 
the 100 torr case using a driver pressure of 3.45 MPa 
(500 psi) and a driver temperature of 330 K. It is 
seen from the one-dimensional simulations that the 
300 K driver temperature would give better agreement 
with experiment but the expensive axisymmetric sim- 
ulations were done before the experimental data was 
gathered. Two different meshes are used for the ax- 
isymmetric investigations both of which contained 700 
points along the length of the tube (350 each in the 
driver and driven sections) and 112 points between 
the tube centerline and an outer wall. The two meshes 
differ in the grid spacing used at the wall. The points 
were exponentially clustered near the wall with a min- 
imum spacing of .075 mm at the wall for one case and 
.0375 mm for the other case. These solutions will be 
referred to as the “coarse” and “fine” solutions, respec- 
tively. Both cases used a larger wall spacing of .15 mm 
in the driver section and in the first half of the driven 
tube. The specified finer spacing was present in the 
final .75 meters of the driven tube after being linearly 
reduced from the larger spacing. This approach eased 
the computational cost by allowing a larger time step 
early in the solution and also by allowing the same 
starting solution to be used for both computations 
since the grid was identical for approximately one half 
of the driven tube. As mentioned before, the solutions 
assumed an initially planar contact discontinuity and 
laminar viscosity. A constant wall temperature of 300 
K was also specified. 

Figure 7 shows computed density and velocity con- 
tours of the end wall region shortly after the reflection 
of the incident shock on the fine grid solution (the con- 
tact discontinuity is still far upstream at this point in 
the simulation). The contours are very similar to the 
schematic diagram in Figure 2. The boundary layer 
ahead of the reflected shock is quite thin and contains 
approximately 16 points. There is a region of sepa- 
rated, recirculating flow under the shock bifurcation. 
A slight curvature of the normal shock which is ob- 
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served in experiments is also perceivable in the simu- 
lation. Along the wall in the stagnation region there 
is a layer of higher velocity fluid that curls up when it 
encounters the end wall. This feature compares quite 
favorably with the schlieren pictures of Kleine ei ai 12 

Figure 8a shows temperature contours of the end 
wall region a short time before the the reflected shock 
interacts with the contact discontinuity. The driver 
gas is easily distinguished with temperature because 
it is much cooler than the driven gas. Note that 
the contact discontinuity is still planar. Because the 
boundary layer in this portion of the tube has had a 
longer time to develop compared to the end wall re- 
gion, it is thicker and contains approximately 30 mesh 
points. The reflected-shock/boundary layer interac- 
tion has significantly changed since the time at which 
Figure 7 was recorded. Although the familiar shock 
bifurcation is still present, there is now a secondary 
shock structure behind the original reflected shock. 
This phenomenon appears similar to the pseudo-shock 
structure mentioned in the introduction. The sepa- 
rated flow region under the foot of the reflected shock 
has expanded so that it extends from the front foot of 
the shock bifurcation to the area under the secondary 
shock. Even with the modified structure of the pseudo- 
shock there is still a wall jet originating at the oblique 
shock of the bifurcated foot and extending all the way 
to the end wall. The jet is visible with temperature 
contours because the temperature rise experienced by 
the gas passing through the oblique shocks of the bifur- 
cated foot is 100-200K less than the rise experienced 
by the gas passing through the normal reflected shock. 

Figure 8b shows temperature contours a short time 
after the reflected shock interacts with the contact dis- 
continuity. The figure indicates that the bifurcated 
shock structure significantly deforms the contact dis- 
continuity. Hot driven gas in the separated region un- 
der the shock bifurcation has been carried into the 
driver gas while the driver gas that passed through the 
oblique shocks of the foot region has begun to pene- 
trate into the driven gas (i.e. it is part of the wall jet). 
There is also the additional feature of a reflected shock 
due to the overtailored conditions. 

Figure 9 contains temperature contours from a time 
later than in Figure 8. The flowfield has become even 
more complicated. The penetration of the driver gas 
along the walls is now clearly seen and the bubble of 
separated flow under the bifurcated shock has carried 
hot driven gas with it far from the end wall. 


Figure 10a, b repeats the experimental data and the 
1-D results with the 330 K driver of Figure 6 and 
adds the results from the two axisymmetric simula- 
tions. Figure 10a is from .445 m upstream of the 
end wall and Figure 10b is .025 m from the end wall 
(the order is reversed compared to Figure 6). As was 
mentioned for Figure 6b, the 1-D formulation is un- 
able to account for the shock/boundary layer interac- 
tion and predicts a slower reflected shock arrival time 
at the pressure transducer location. The axisymmet- 
ric solutions include this phenomenon (with a laminar 
shock/boundary layer interaction) and predict a faster 
reflected shock speed and closer agreement to the ex- 
periment. Strehlow and Cohen 17 also reported that 
the reflected shock speed was increased by the pres- 
ence of the shock bifurcation. By considering the rela- 
tive change produced in Figure 6b by using a driver at 
300 K rather than at 330 K, it is presumed that even 
better agreement with experiment would result by us- 
ing the lower driver temperature in the axisymmet- 
ric simulation. Another factor that may influence the 
comparison between the experimental and simulated 
reflected shock speeds is the character of the bound- 
ary layer. The presence of a turbulent boundary layer 
in the experiment would cause a growth of the bifur- 
cated shock different than the growth computed with 
the laminar viscosity assumption. 

The general features of the pressure traces obtained 
near the end wall from the one-dimensional and ax- 
isymmetric solutions in Figure 10b are fairly similar. 
One difference is the presence of small oscillations in 
the computed axisymmetric pressure traces (particu- 
larly after the tailoring wave arrives). When looking 
at the many flow features created in the stagnation 
region in Figures 8 and 9, it is easy to understand 
that any reflected waves interacting with these fea- 
tures would produce very complicated pressure pat- 
terns. Figure 10b also allows an opportunity to com- 
pare the experimental and computed pressure traces 
“created by the passage of the bifurcated shock. An 
enlargement of this region is presented in Figure 11. 

Figure 11 reveals that there are several significant 
differences between the experimental and the compu- 
tations. The largest discrepancy is the large computed 
pressure drop under the bifurcated foot where the ex- 
periments show more of a pressure plateau. This dis- 
crepancy is reinforced by the fact that going from the 
coarse grid to the fine grid increased the size of the 
computed pressure drop. The explanation for this dif- 
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ference is not known. It is interesting to note, how- 
ever, that the experimental pressure traces presented 
by Sanderson 15 show a drop in pressure under the sep- 
arated region (comparable to the current simulations) 
and that similar behavior is sometimes observed in the 
current experimental apparatus under different initial 
conditions. These observations support the physical 
nature of a pressure drop as well as the the ability of 
the transducers to capture such a phenomenon. An- 
other difference between the experiment and the ax- 
isymmetric simulations is the size of the pressure over- 
shoot after the separation bubble. Part of the dis- 
agreement can be attributed to the width of the pres- 
sure transducer (.21 inches). Both the pressure drop 
mentioned previously and the overshoot are created by 
features which are smaller than the width of the trans- 
ducer so that some smearing of the experimental data 
is probably occurring. To illustrate possible influence 
of the transducer size, Figure 12 repeats the pressure 
trace in Figure 11 for the fine grid axisymmetric so- 
lution and also contains a pressure trace created by 
integrating the computed pressure over the surface of 
the transducer and dividing by the total transducer 
surface area. Some smearing in the computed traces is 
produced by this averaging but significant differences 
between the computations and the experiment remain. 

There is room for improvement in the current sim- 
ulations. While sufficient resolution has been achieved 
in the boundary layers for computing laminar velocity 
profiles, the thermal boundary layer is not sufficiently 
resolved. Future work will explore more formally ex- 
amine the time accuracy of treating the viscous terms 
implicitly. The accuracy of treating inviscid terms im- 
plicitly will also be investigated. It is hoped that if 
an implicit time step at a CFL number greater than 1 
based on the inviscid terms can be used for the small 
cells near the wall (while still retaining time accuracy 
in the regions away from the wall), then a much finer 
wall spacings and lower computational cost may be 
possible. Better resolution may explain some of the 
differences seen in this work between the experiments 
and computations. 

The knowledge gained through the present research 
will be used to design experiments so that data can be 
acquired for validation of the predicted pseudo-shock 
structures and the arrival of driver gas in the stagna- 
tion region. Initial attempts at measuring this latter 
phenomenon will be done with heat-transfer measure- 
ments. Further studies will also allow many of the 


phenomena investigated analytically and experimen- 
tally by earlier works to be investigated numerically. 
These topics include shock bifurcation growth rates 
and the effects of turbulent boundary layers. 

Conclusions 

Axisymmetric simulations of the NASA Ames 
electric-arc driven shock tube have been done which in- 
clude the wall boundary layer. The computed results 
have been compared with experimental data. These 
simulations have allowed the wall jet created by the 
reflected-shock/boundary layer interaction to be inves- 
tigated numerically for the first time. These simula- 
tions support earlier analytical and experimental work 
which indicate that this mechanism can contribute to 
the reduction of the usable test time by allowing the 
driver gas to contaminate the stagnation region. It was 
also shown that even before any driver gas contamina- 
tion, the wall jet creates a layer of lower temperature 
gas on the shock tube walls which might effect the in- 
terpretation of temperature measurements in the stag- 
nation region. These results reinforce the fact that 
the stagnation region in shock tubes and shock tun- 
nels is a region of complex flow phenomena. There 
were discrepancies between the experimental and com- 
puted pressure traces under the bifurcated foot of the 
reflected shock. The differences may be due to a lack 
of grid resolution or turbulent boundary layers. 

The present results show that multi-dimensional 
simulations of shock tubes and tunnels are expensive 
but also valuable. An implicit treatment of the viscous 
terms was required to make the simulations practical. 
The lessons learned from this initial work are being 
used to determine future experimental and computa- 
tional directions for the purpose of better understand- 
ing shock tube flows. Future experimental work will 
include gathering data at more locations to help ver- 
ify the prediction of pseudo-shocks and the use of heat 
transfer and visual data to further study the driver gas 
contamination process. 
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Fig. 2 Schematic diagram of the interaction of the Fig. 3 Schematic diagram of the pressure trace under 
reflected shock with the boundary Layer. a bifurcated shock. 


.0254 m 



End 

wan 


Fig. 4 Geometry and instrumentation for the NASA electric-arc driven shock tube. 
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Fig. 5 Computed x-t diagram of temperature contours for a one-dimensional simulation of the NASA Ames 
electric-arc driven shock tube. 




Fig. 6a Experimental and computed pressure traces Fig. 6b Experimental and computed pressure traces 

.0254 m from the driven tube end wall. .445 m from the driven tube end wall. 
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Density Contours Axial Velocity Contours 



Fig. 7 Density and axial velocity contours showing the computed reflected shock/boundary layer interaction 
(t = 2.74 msec). Compare to the schematic diagram in Figure 2. 




Fig. 8a, b Computed temperature contours before and after the interaction of the bifurcated reflected shock 
wave with the contact discontinuity: a) t = 3.00 msec b) 3.06 msec. 
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Fig. 9 Computed temperature contours after the interaction of the reflected shock with the contact discontinuity 
(t = 3.25 msec). 




Fig. 10a Experimental and computed pressure traces Fig. 11 Enlargement of Figure 10b showing the pas- 

.445 m from the driven tube end wall. sage of the bifurcated reflected shock. 



Fig. 10b Experimental and computed pressure traces Fig. 12 Computed pressure trace showing the effect 

.0254 m from the driven tube end wall. of the pressure transducer surface. 
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Abstract 


The flowfield around a proposed lunar return aerobrake is computationally examined 
using both axisymmetric and three-dimensional formulations. Aerobrake configurations with four 
panel-shaped braking surfaces at 50°, 55°, and 60° to the freestream are calculated and the 
complex nature of the three dimensional flowfield is examined. It is shown that, despite the 
complexity of the three dimensional flowfield, the axisymmetric solutions give results very similar 
to the three-dimensional solutions for total drag coefficient. The full paper will examine the effect 
of changing the width of the braking surfaces on the relationship between the axisymmetric and 
three-dimensional calculations for drag. 


Introduction 

Future space exploration holds many exciting possibilities not the least of which are 
manned and unmanned missions to other planets. For this type of exploration, a lunar base may 
be employed as a way station for personel and equipment between the target planet and Earth. 
Transportation from a lunar base to Earth can be accomplished by an aerobraking vehicle. The 
lunar return aerobrake concept examined here is essentially a blunt nosed cylinder that uses four 
large, flat panels as braking and control surfaces, Fig. 1 . The nose of the aerobrake is an ellipsoid 
so that drag from the nose region is maximized and heating is minimized. The radius of the body 
increases linearly from the beginning of the main body .where the nose ends, to the end where 
the panels attach. The main body is flattened where the panels attach such that the width of the 
main body in that region is equal to the nose diameter. The tapering of the lower portion of the 
panels and the flattening of the main body allow for a long hinge line so that the bending loads on 
the panels can be distributed along the hinge without unduly reducing the volume of the vehicle. 
The four panels are extended at angles of attack to the flow and are independently adjustable so 
that desired angles of attack and lift to drag ratios can be achieved (L/D = 0.3 for lunar return). The 
panels are large so that the ballistic coefficient (M/CdA) is reduced thereby reducing the total 
heating to the vehicle. However the width of the panels is limited so that they can be stowed at 0° 
deflection while inside the launch vehicle. The payload capacity of the aerobrake is estimated to 
be 10,000 kg and the total volume to be 150 m3. The vehicle is designed such that it could be 
launched from Earth aboard the proposed National Launch System (NLS), and it can also fit within 
the Space Shuttle payload bay (4.57 m x 18.28 m). 

A large portion of the aerobrake trajectory will be dominated by real-gas effects. These 
effects can be approximated for initial design purposes by using an ideal-gas assumption which 
incorporates an effective value of specific heat ratio, y. This method has been shown to be 
suitable for the calculation of drag and moment coefficients of bodies in real-gas flows. 1 > 2 The 
trajectory point chosen for the study corresponds to the peak pressure trajectory point, Table 1. 

At these freestream conditions the stagnation pressure is nearly 0.4 atm. and the flow is assumed 
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to be in chemical and thermal equilibrium. Therefore, the effective value of y, y=\2\, is chosen to 
be the specific heat ratio that corresponds to equilibrium conditions at the given speed and 
altitude. 


Numerical Methods 

The laminar flow about an aerobrake with all four panels at three different orientations, 
consisting of 50°, 55°, and 60°, is examined computationally. Because of the flat panels and the 
flattened portion of the main body, the aerobrake has only one octant symmetry (45° from the 
middle of a panel to the middle of a gap region). Therefore, a code that is capable of solving a 
three-dimensional (3-D) problem is required. However, the geometry is very nearly axisymmetric 
and can be estimated by creating a body of revolution using a line that runs from the nose of the 
aerobrake down the center of a panel. For this approximated geometry, an axisymmetric solver 
can be used. Solutions are computed for both the 3-D geometry as well as the axisymmetric 
geometry. Although the axisymmetric case runs much faster it has the disadvantage that the gaps 
that exist between the panels are neglected and the braking surface becomes a solid cone 
shape. Therefore, the effect of the gaps on the aerobrake flow cannot be examined using the 
axisymmetric geometry. 

Two different CFD codes were used to calculate the axisymmetric flow around the 
aerobrake. The first code is an axisymmetric solver developed by Candler.3 This is a finite-volume 
code that uses flux-vector splitting and Gauss-Seidel line relaxation to solve the Navier-Stokes 
equations. This code has been used previously for similar studies of hypersonic drag brakes. 2 
The second code used is a 3-D solver, FL3D, developed by Venkatapathy.4 This is a finite 
difference code that uses a three-dimensional, implicit time marching Navier-Stokes solution 
procedure. This code has previously been applied to axisymmetric and three-dimensional nozzle- 
plume and blunt body flows. The numerical method is a LDU-ADI scheme with Roe's averaging 
and MUSCL differencing. High order spatial accuracy is achieved by constructing MUSCL 
extrapolation of flow variables with a differentiable limiter. Roe averaged fluxes are constructed 
and the flux-difference from the higher order fluxes are used in this formulation. The LDU-ADI 
algorithm is a diagonal algorithm requiring minimal CPU per iteration and is applicable to steady 
flows, and the ADI sweeps of this formulation are quite efficient. Both the codes are ideal-gas 
solvers which can incorporate an effective value of y. 

Since no experimental data for this type of flowfield exist, the two codes were compared 
against each other to assure that the flowfields are calculated in a consistent manner. The codes 
were used to calculate the flowfield around an axisymmetric aerobrake with panels at 50°. The 
general flow structure, including the size of the separation bubble caused by the panels, 
compared well, and the differences in calculated surface pressures were less than 1%. The size 
of the separation region is important because it impacts the efficiency of the braking panels. The 
pressure within a separated region remains essentially constant even when part of that region is a 
compression surface such as the base of a panel. Therefore, the pressure increase that would 
exist if no separation occurs is not realized due to the separation. 

For each set of panel angles, an axisymmetric and a 3-D solution were calculated. 

Although the 3-D solver could be used to calculate the axisymmetric solutions as well as the 3-D 
solutions, the axisymmetric solver was used because of its speed and convenience. Each 
axisymmetric solution was started from the same initial conditions and required approximately IQ- 
15 minutes of CPU time on a Cray Y-MP. The first 3-D solution, the 50° panel case, was started 
from a blunt cone solution and required approximately 1.75 hours of CPU time. Subsequent 3-D 
solutions were started from the most similar preexisting 3-D solution. For example, the solution 
from the 50° panel case was used as the initial conditions for the 55° panel case. This reduced 
the CPU time to approximately 45 minutes. 

The grids used are 61 x47x63 for the 3-D calculations and 61 x63 for the axisymmetric 
calculations. A representative axisymmetric grid is shown in Fig. 2. For the axisymmetric 
calculations, a solution-adapted 5 grid technique was used to refine important flow features. This 
technique will also be used for future 3-D calculations. For the 3-D calculations, symmetry is taken 



advantage of and only one octant (45°) of the full 3-D geometry is computed. The surface grid for 
the 3-D calculations is shown in Rg. 3. Because the nose of the aerobrake is truly axisymmetric at 
the zero angle of attack considered here and is not influenced by the flow downstream, there is 
no need to compute the nose region using a 3-D formulation. Therefore, the portion of the 
flowfield that is actually computed by the 3-D solver begins at the main body. The inflow of the 3-D 
calculation is fixed to the values from the axisymmetric solution for the nose region. This reduced 
the grid size for the 3-D calculations to 46x47x63. 

The assumption is made that the base flow region has very little affect on the total drag on 
the vehicle, since the pressure in that region is relatively low. However, because of the gaps, the 
base flow may have some effect on the main body. It is not clear what these effects may be, and a 
study of the full flowfield including the base flow is reserved for future work. 

Results and Analysis 

A comparison of the axisymmetric and 3-D calculations of surface pressure is made for an 
aerobrake with panels at 50°, Fig. 4. Although the axisymmetric solution shows a somewhat 
higher peak pressure on the panel, the surface pressure profiles look very similar indicating that, 
at least along that one line of the 3-D geometry, the axisymmetric solution comes very close to the 
3-D solution. However, the limitations of the axisymmetric solution become apparent when 
simulated surface oil flow for the 3-D solution, Fig. 5, is examined. The surface pressure on the 
body in Fig. 5 is indicated by colors; dark blue is the lowest pressure, essentially freestream, and 
white the highest. The oil flow is show in yellow. The separation on the main body is largest in the 
center of the panel regions and decreases only slightly to the smallest size in the centers of the 
gap regions. The patterns of the oil flow indicate that a complex, multiple vortex structure exists 
within the separated region. This structure will be examined to a greater extent in the full paper. 

Particle traces constrained to the plane that lies on the center of a panel are shown in 
yellow in Fig. 6. The particle traces indicate the size of the separated region and show more 
clearly than the surface oil flow the extent of the separation on the panels. Almost half of the 
panel length is contained within the separation. The effect of the separation on the panels is 
demonstrated by the painted surface pressure contours which show a very low pressure within 
the entire separated region. A cut away of the outer shock region is shown in this figure as a solid 
surface. The effect of the separation on the outer shock surface is indicated by the slight bulging 
of the surface in the region of the separation. Because of the flare-gap geometry, the shock 
structure in the flare region is truly three dimensional. It retains a conical shape in the gap regions 
and bulges out in the panel regions where it reflects off the panels. The shock stand off distance 
on the panels is greatest at the center of the panels. As the flow expands around the edge of the 
panels, the shock curves and the stand off distance decreases. This causes the peak pressure 
on the panels to shift off the center line to a position approximately midway between the panel 
center and the edge. 

Particle traces constrained to the plane that lies in the center of a gap region are shown in 
Fig. 7. The particle traces terminate at the end of the computational domain. The separation in 
the gap region remains relatively large and, as is indicated by the traces, the recirculation region 
extents into the base area where the flowfield is not calculated. The flow at the outflow of the 
computational domain remains supersonic even within the separation, therefore, the 
extrapolation boundary conditions applied at the outflow region are valid. However, the effect of 
the base flow region on the main body is not captured by this representation. 

The area of a circle with a radius equal to the outer radius of the base of the aerobrake 
main body is used as the reference area for calculating drag coefficients For the axisymmetric 
solutions, the total Cd is computed by calculating the drag assuming the panel to be solid and 
then subtraction out the portion of the drag that is produced where the gaps occur in the 3-D 
geometry. Since the 3-D flares are flat and not rounded as is assumed for the axisymmetric 
solution, the flare size for the axisymmetric solution is determined by assuming that the width of 
the flare is the base of a triangle with its vertex at the symmetry line. The angle opposite the base 



of the triangle is then used to determine the arc length of the "axisymmetric flare". The tapered 
portion of the 3-D flare is assumed to be straight for the axisymmetric calculations. 

A comparison of the total Cp calculated for the axisymmetric and 3-D solutions for 50°, 
55°, and 60° flare angles is shown in Fig. 8. The most notable feature of this comparison is that, 
although the axisymmetric solution is incapable of computing the complex flowfield revealed by 
the 3-D solutions, the total drag calculated for the axisymmetric and 3-D solutions are very similar. 
This suggest that, for this type of geometry, it is not necessary to compute a full 3-D solution in 
order to calculate drag. Presently, calculations are being performed to compare different "f lare-to- 
gap" ratios in order to determine whether the similarity between the axisymmetric and 3-D drag 
calculations is maintained. For these calculations, the flat-flare 3-D geometry is replaced by a 
geometry that would be axisymmetric if there were no gaps between the braking surfaces. This is 
done so that a one to one comparison of drag can be made between the 3-D and axisymmetric 
solutions. The results form these calculations will be given in the full paper . 

Conclusions 

The 3 D calculations of a proposed lunar return aerobrake show that a large portion of the 
main body as well as the braking surfaces are contained within a large, complex separated region. 
Despite the complexity of the 3-D flowfield, axisymmetric solutions predicted total drag 
coefficients very close to those calculated from the 3-D solutions for aerobrakes with panels at 
50°, 55° and 60°. It is presently being determined whether different "flare-to-gap" ratios, i.e. 
different flare widths, will change the relationship between the axisymmetric and 3-D calculations 
of drag. 
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altitude 

56.6 Km 

velocity 

9082 m/s 

density 

4.751x10-4 

pressure 

35.0 N/m 2 

temperature 

256.6 K 

wall 

temperature 

2700 K 


Table 1 Freestream Conditions 
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Numerical Study of the Transient Flow in the Driven Tube and the 
Nozzle Section of a Shock Tunnel 

Susan Tokarcik and Jean-Luc Cambier 
Eloret Institute 

ABSTRACT 

The initial flow in the NASA Ames 16" shock tunnel is examined numerically. 

A finite-volume total variation diminishing (TVD) scheme is used to calculate 
the transient flow in the region of the driven-tube/nozzle interface. The flow 
is assumed to be inviscid and non-reacting. Two types of calculations are 
presented: a series assuming that the flow is axisymmetric, and a calculation 

assuming the flow is two dimensional. The differences between the 2-D and 
the axisymmetric flows are examined. The axisymmetric solutions predict that 

a vortical structure is formed in the driven-tube/nozzle interface region and 
propagates into the driven tube along with the reflected shock. Several 
parameters are examined in order to determine their effect on the vortex 
structure: 1) varying throat size, 2) changing the shape of the juncture 

between the driven tube and the nozzle, and 3) changing the shape of the 
driven tube's end wall. 

INTRODUCTION 

The Ames 16" facility, shown schematically in Figure 1, can be considered as a 
typical example of a shock tunnel for hypersonic flows. The shock tunnel is 
about 70 meters long and is composed of a driver section (17" diameter) filled 
with a light combustible mixture, a long driven tube (filled with the test gas at 
low pressure), a supersonic nozzle, and finally a test section. The test times of 
this facility are typically 5 to 20 milliseconds, and the time for the main shock 
to propagate from the main diaphragm to the end of the driven tube is on the 
order of 7 millisec. The shock partially reflects at the driven tube end-wall 
bursting the second diaphragm and initiating flow in the nozzle. The reflected 

shock eventually interacts with the contact discontinuity (CD) which separates 
the test gas from the driver gas. More detailed information about shock tunnel 
operations can be found in ref. [1]. The test time is measured from the time the 
flow conditions in the nozzle stagnation region become steady, until the CD 
reaches the end of the driven tube, and the nozzle flow becomes contaminated 
with the driver gas. 

The initial, unsteady flow in the shock tunnel is important because it directly 
affects the quality of the resulting test flow. Because the initial flow is 
difficult to examine experimentally, numerical analysis is used to simulate the 
start-up process. Recent numerical studies, ref. [2] and ref. [3], have predicted 
that the flow in the region of the driven tube and nozzle interface is quite 
complex. A particularly important feature is a vortex system which is formed 
in the center of the driven tube. This vortex may have an impact on the 
quality of the test flow as well as the overall test time; it is therefore examined 
in more detail in this study. 

NUMERICAL ANALYSIS 

A Finite-volume, explicit, total variation diminishing (TVD) scheme developed 
by Cambier is used to calculate the flow which is assumed to be inviscid and 


non-reacting. The code has first order as well as second order spatial accuracy 
capabilities. Most calculations were performed using second-order accuracy; 
calculation which used first-order accuracy are noted as such. This code has 
been employed previously to calculate similar transient flows, see ref. [2], 

90 deg Expansion Corner Flow 

A validating exercise was performed to show that the code is capable of 
adequately resolving transient flows which are highly vortical in nature. A 
schematic sketch of the shock tube experiment, which was conducted by 
Mandella [4], and Mandella and Bershader [5] is shown in Figure 2. The shock 
tube, which has a 5cm x 5cm cross section, produced a planar shock wave with 
a Mach number M s =2.0. The end of the shock tube was open to ambient 
conditions. Two parallel plexi-glass plates were attached to the end of the 
shock tube in order to keep and shock wave and the rest of the flow field two 
dimensional. The planar shock wave starts to diffract as soon as it leaves the 
shock tube and forms a curved shock wave. The experimental interferogram 
—(Figure 3) show that, as the shock wave propagates into the test section, 
vortices are formed in the vicinity of each of the 90 deg expansion comers. 
Also seen are the contact surface, slipstream, and Prandtl-Meyer fan which 
interact with the vortices. 

Figure 4 shows computed density contours at approximately the same time as 
the experimental interferogram. Comparing the two figures demonstrates that 
the computations are able to predict all the major flow features mentioned 
above with good accuracy. 

Shock Tunnel Driver/Nozzle Flow 

The starting conditions used for all the shock tunnel calculations correspond 
to the point in time when the main shock has reached the driven tube's end 
wall and the secondary diaphragm has ruptured. The flow conditions in the 
driven tube at this initial point are calculated taking into account chemistry 
and correspond to the conditions for a main shock with a speed of 12 km/sec 
and a driver tube pressure of 6000 lbs/sq. in: P=50.76 ATM, T= 3674 K, U=2710 
m/s, N=0.6042xl0' 4 , 0=0.04091, N2=0.7341, 02=0.1456, N0=0.07934 (unit of mole 
fraction). The initial conditions in the nozzle correspond to air with a 79% N 2 
and 21% O 2 mixture at a pressure of 1.3xl0‘ 4 ATM and a temperature of 298 K. 

As mentioned earlier, chemical reactions are not considered in this study, and 
the flow is assumed to be chemically frozen. 

As was mentioned earlier, the initial flow in the region of the driven- 
tube/nozzle interface is quite complex. Numerical calculations have predicted 
that a vortex forms at the center line of the driven tube and propagates with 
the reflected shock into the driven tube. The vortex can cause the reflected 
shock to bulge at the center line (Figure 5) which may cause premature 
contamination of the test gas by the driver gas thus decreasing the test time. 

Several aspects of the calculated transient flow are examined: 1) the effects of 

the numerical method on the resulting flow field, 2) the effects of geometry, 3) 
the effects of time (how the flow field changes as the reflected shock and 
vortex structure continues to propagate into the driven tube and 4) an 
axisymmetric versus two-dimensional flow field. 


(Note: the full paper will also examine the interaction of a bulging, reflected 
shock with a CD.) 

Numerical Effects 

Two components of the numerical technique are examined: 1) grid effects and 

2) first order versus 2nd order spatial accuracy. Figure 6a shows a second- 
order solution on a uniform grid at a time t=80 p.sec after shock arrival at the 
driven lube's end wall (80 psec from starting the calculation). This figure 
shows that the reflected shock has formed a coned shaped at the centerline. 
Figure 6b shows a solution at the same time (and same accuracy), however the 
aspect ration of the grid is AX:AY=3:1. In other words, the grid is clustered at 
the axis. Comparing these two figures it becomes evident that the jetting 
effect at the axis has been aggravated by the aspect ratio of the grid. This 
peculiar formation can also be seen in the results of P. Jacobs, ref. [3]. This 
structure is possible if an intense and high velocity jet of gas is produced and 
maintained on the axis, which is a highly singular and unphysical behavior. 

_ One explanation of this phenomenon, which is described in ref. [2], is that the 
jetting is caused by the numerical representation of the singular axis in the 
axisymmetric formulation. It was assumed in ref. [2] that the problem arose 
from the axisymmetric pressure correction term which is not part of the 
monotonic (TVD) fluxes and acts as a non-conservative momentum source that 
behaves as l/(Ar) near the axis. Therefore clustering points near the axis 
(thus accentuation the l/(Ar) term) aggravates the jetting. However, further 
study has shown that the jetting phenomenon can also occur off the axis as^ 
well as in solutions which have two-dimensional formulations. Therefore, ~ 
although the pressure correction term is still believed to be a contributing 
factor in the axis problem, it has become clear that there is still an, as yet, 
unknown cause for the jetting in other cases. 

Despite not fully understanding the cause of the problem, it has been 
determined that coarsening the grid in the radial direction alleviates the 
problem in all instances. Figure 6c show the same case as computed in Figures 
6a and 6b, however the grid aspect ratio has been changed so that AX:AY=3:1 
(near the axis). Comparing Figures 6a, 6b and 6c shows that the unphysical 
jetting effect, which is present in both the AX:AY=3:1 and the uniform grid, is 
eliminated in the third solution by the proper choice of grid aspect ratio. 

Another way to somewhat alleviate the jetting problem is to compute the 
solution using only first order spatial accuracy. This essentially causes the 
numerical scheme to become more dissipative which in turn mitigates the 
jetting effect. Using a first order formulation does have some serious 
consequences, however. Figures 7a and 7b compare two solutions, each at 
t=200 psec and each on the same AX:AY=1:3 grid. The first solution employed 
first-order spatial accuracy while the second used second-order accuracy. 
These figures show that as the reflected shock continues to propagate down 
the driven tube, the dissipative, first-order scheme loses any trace of the 
vortex structure while the second-order solution clearly shows that the vortex 
structure not only remains but has increased in size. Therefore, it is clear 
that the study of this transient flow must be conducted using at least second- 
order accuracy and employing grids which are coarsened in the radial 
direction. All flow solutions presented in the remainder of the study employ 
the latter two criteria. 


Effects of Geometry 


Assuming that all major numerical effects have been removed from the 
problem, the effects of changing the geometry in the driven-tube/nozzle 
interface region is now studied. Several parameters are examined: 1) 
changing the throat size, 2) changing the shape of the juncture between the 

two sections, and 3) changing the shape of the driven tube's end wall. 

The area ratios of the nozzle throats examined are A/A*=100, A/A*=190, and 
A/A*=570. Figures 8a, 8b, and 8c show temperature contours for these three 
throat sizes at t=120 nsec. A comparison of the temperature contours indicates 
that, while the flows within the nozzle sections are somewhat different in each 
case (in particular, the shock interaction in the subsonic portion of the throat 
is intensified by the smaller throats), the vortex structures in the driven tubes 

are identical. Therefore, it is concluded that the size of the throat has at most 

a minimal affect on the formation of the vortex in the driven tube. 

The affect of the size of the aperture at the end of the driven tube is now 
examined. The size of the hole in the driven tube was increased by 35%, Figure 
9, and the flow was recomputed for the A/A*=190 throat. The temperature 
contours show that for this case the flow in the nozzle section is unaffected, 
and the flow is the driven tube is somewhat affected. The central vortex 
structure is unchanged, however the length of the "arms" on either side of the 
central vortex have increased causing an overall size increase of 16% over the 
case discussed above (which is significantly less than the 35% increased in 
aperture height). From this it can be concluded that increasing the aperture 
height increases the overall size of the affected area of the reflected shock, 
however the increase is not one to one. 

Another change that can be made to the driven-tube/nozzle juncture is 
rounding the interface. Rounding the interface resulted in an aperture 

which was slightly larger that the 35% increase case discussed above. Again, 
the overall structure and size of the central vortex was unchanged, however 
the size of the affected area increased significantly (43% larger than for the 
original geometry). Another important difference was that the center vortex 
observed in this case protruded slightly further into the driven tube than for 
any of the other cases (all of which protruded to the same distance). Also, the 
flow field in the nozzle section is affected in this case. In the subsonic portion, 
the shocks emanating from the juncture comers are no longer present, and 
compression waves have appeared downstream of the throat. 

The last geometric change that is examined is altering the shape of the end- 
wall of the driven tube. Figure 10 shows how the end-wall has been angled 
outward from the juncture so that it is no longer vertical and, again, the 
juncture between the two section is rounded. This flow field (t= 120 p.sec) looks 
considerably different than those examined so far; there are strong shocks in 
both the subsonic and supersonic portions of the nozzle, and the character of 
the vortex structure has changed. The overall vortex structure is much larger 
than for the original geometry, and the central structure extends 
significantly further into the driven tube than for any of the other cases. 

Also, the central vortex itself is significantly larger than in the other cases 
where the height of the overall structure changed but the central structure 
remained virtually unchanged. 


For the bulging of the reflected shock to be of importance it would have to 
persist as the shock continued to travel up the driven tube to the CD where the 
bulge would cause the driver gas to be prematurely mixed with the test gas. 

The original A/A*=190 case was continued for a significant time (t= 400psec), 
and it was found that the vortices did appear to dissipate as the reflected shock 
continued to propagate down the driven tube. Therefore, the likelihood of 
premature mixing will depend on how much time elapses before the CD and the 
reflected shock interact. 

A question that arises is whether the formation of the vortex is an 
axisymmetric phenomenon. Figures 11a and lib show temperature contours 
for a 2-D calculation using the original A/A*=190 geometry at t=120 psec and 
t= 1 60 psec. For the two-dimensional cases, the centerline vortex does not form; 
furthermore, the portion of the reflected shock near the centerline actually 
lags behind the rest of the shock. Even at t=160 psec the centerline region of 
the shock has not caught up with the rest of the shock. 

The behavior of the flow field for a three-dimensional formulation may be 
presented in the full paper. 


CONCLUSIONS 

It was found that refining the grid in the radial direction (such that Ax :Ay > 1 
produced an unphysical, jetting behavior in the solution. A partial 
explanation was given, and changing the grid aspect ratio such that Ax: Ay > 3 
was presented as a solution to the jetting problem. Spatially first order 
accurate solutions were found to be too dissipative to resolve the vortex 
structure in the driven tube. 

Several geometry changes were made to both the nozzle section and the driven 
tube to determine their effect on the flow field. It was found that 1) changing 
the throat area had an unobservable effect on the flow in the driven tube and 
relatively minor effects on the nozzle flow, 2) increasing the height of the 

aperture in the driven tube increased the size of the overall vortex structure; 

however the size increase was not one to one and the size of the core structure 
was unchanged, 3) rounding the juncture between the nozzle and driven 

sections caused the size of the vortex structure to increase markedly, and also 
caused the reflected shock to bulge slightly more 4) changing the shape of the 
driver tube's end-wall such that it was no longer vertical caused a more 
intense vortex to form in the driver tube which caused a significantly larger 
bulge in the reflected shock, and 5) the vortex structure eventually dissipated 
at the reflected shock continued to propagate into the driven tube. 

Solutions assuming that the flow was two dimensional showed no that no 
vortex existed in the driven tube, and the reflected shock near the centerline 
actually lagged behind the rest of the shock instead of protruding ahead of it. 
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Particle Traces Colored by Pressure 
Demonstrate the Existence of a Vortex 
System At the Center of the Driven 
Tube which Causes the Center Portion 
of the Reflected Shock to Bulge 
Upstream. 



Temperature Contours at 
time=80psec Computed 
on a Uniform Grid with 2nd 
Order Spacial Accuracy 



Fig. 6b Temperature Contours at 
time=80psec Computed 
on a Grid with Aspect Ratio 
Ax:Ay=3:l 


2nd Order 
X:Y=3:1 



Fig. 6c Temperature Contours at 
time=80|isec Computed 
on a Grid with Aspect Ratio 
Ax:Ay= 1 :3 




Fig. 7a Temperature Contours at 
time= 1 20(isec Computed 
on a Grid with Aspect Ratio 
Ax:Ay=3:l and 1st Order 
Spacial Accuracy 


Fig. 7b Temperature Contours at 
time= 1 20(isec Computed 
on a Grid with Aspect Ratio 
Ax:Ay=3:l and 2nd Order 
Spacial Accuracy 


2nd Order 
X:Y=1 :3 



Fig. 8a Temperature Contours at 
time= 1 20psec Computed 
on a Grid 'with Aspect Ratio 
Ax:Ay=3: 1 and 2nd Order 
Spacial Accuracy for a Large 
Throat Size 


Fig. 8b Temperature Contours at 
time= 1 20|isec Computed 
on a Grid with Aspect Ratio 
Ax:Ay=3:l and 2nd Order 
Spacial Accuracy for a Medium 
Throat Size 


Fig. 8c Temperature Contours at 
time= 1 20psec Computed 
on a Grid with Aspect Ratio 
Ax:Ay=3:l and 2nd Order 
Spacial Accuracy for a Small 
Throat Size 



Fig. 9 Temperature Contours at 
time= 1 20psec Computed 
for a Geometry with a 35% 
Larger Aperature in 
the Driven Tube End Wall 



Fig. 10 Temperature Contours at 
time= 1 20psec Computed 
for a Driven Tube with a Non- 
Vertical End Wall 



Fig. 11a Temperature Contours at 
time= 1 20|isec Assuming 
That the Flow is Two 
Dimensional 



Fig. 1 lb Temperature Contours at 
time=l 60psec Assuming 
That the Flow is Two 
Dimensional 
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Numerical simulations are used to investigate periodic combustion instabilities observed in ballistic-range 
experiments of blunt bodies flying at supersonic speeds through hydrogen-air mixtures. The computations are 
validated by comparing experimental shadowgraphs with shadowgraphs created from the computed flowfields 
and by comparing the experimentally measured instability frequencies with computed frequencies. The numerical 
simulations use a logarithmic transformation of the species conservation equations as a way to reduce the grid 
requirements for computing sbock-induced combustion. The transformation is applied to the Euler equations 
coupled to a detailed hydrogen-air chemical reaction mechanism with 13 species and 33 reactions. The resulting 
differentia] equations are solved using a finite volume formulation and a two-step predictor-corrector scheme to 
advance the solution in time. Results are presented and compared for both a flux-vector splitting scheme and an 
upwind TVD scheme. The computations add insight to the physical processes observed in the experiments and 
the numerical methods needed to simulate them. The usefulness of the ballistic-range experiments for the 
validation of numerical techniques and chemical kinetic models is also demonstrated. 


Introduction 

B ALLISTIC-RANGE experiments performed in the 1960s 
and early 1970s 1 " 7 provide excellent data for studying the 
coupling between supersonic fluid dynamics and nonequi- 
librium chemical kinetics as well as for evaluating combustion 
flow codes. In these experiments, small projectiles were fired 
at supersonic speeds into a variety of premixed combustible 
mixtures. Shadowgraphs of the flowfields exhibit two distin- 
guishing features: one is the bow shock ahead of the projectile 
and the other is an energy-release front created by the ignition 
of the heated mixture behind the bow shock. Both features can 
be seen in Fig. 1 in the shadowgraph by Lehr 5 of a steady 
ballistic-range experiment at Mach 6.46. The region between 
the bow shock and the energy-release front is called the induc- 
tion zone, and it exists because there is an ignition delay caused 
by chemical nonequilibrium. The induction zone is character- 
ized by near-constant values (the postshock values) of temper- 
ature, density, pressure, and velocity. The size of the induction 
zone is determined by the fluid speed downstream of the shock 
and the ignition time corresponding to the postshock condi- 
tions. When ignition occurs, the energy is released over an 
interval much shorter than the ignition delay time and appears 
as a discontinuity that is referred to as the energy-release front. 
Across the energy-release front the pressure is nearly constant, 
the temperature rises, and the density drops. This can be seen 
in the shadowgraph where the shift from light to dark across 
the bow shock corresponds to the density increase across the 
shock whereas the shift from dark to light across the energy-re- 
lease front indicates a drop in density. 

Depending on the conditions of the experiment, one will 
observe either steady or unsteady flow. Projectile speeds 
greater than the detonation wave speed tend to induce steady 
flows whereas speeds less than the detonation wave speed can 
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produce unsteady flows. The unsteady flows are characterized 
by two different regimes. One is called the regular regime, and 
the other is called the large disturbance or irregular regime. All 
of the unsteady simulations presented in this work are in the 
regular regime (see Alpert and Toong 7 for more on the large 
disturbance regime). Figures 2a and 2b show an unsteady bal- 
listic-range experiment with the same freestream conditions as 
in Fig. 1 except that the projectile speed is Mach 4.79. The 
shadowgraphs reveal remarkable high-frequency oscillations. 
The frequency of the instability is approximately 720 kHz as 
deduced from the projectile speed and counting the number of 
oscillations that occur over a known distance. Figures 2a and 
2b are from the same ballistic-range shot but show different 
view angles. The view axis of Fig. 2a is perpendicular to the 



Fig. 1 Shadowgraph of a spherical nose projectile moving at Mach 
6.4 6 in a stoichiometric hydrogen-air mixture (courtesy of H. F. Lehr). 
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Fig. 2 Shadowgraph of a spherical nose projectile moving at Mach 
4.79 in a stoichiometric hydrogen-air mixture: a) side view, b) off-axis 
view (courtesy of H. F. Lehr). 


flight axis, and Fig. 2b is off axis and reveals some of the 
three-dimensional structure of the flowfield. Although com- 
plex in many ways, the physics of these ballistic-range flows 
are predominantly driven by reaction kinetics and convection 
phenomena; therefore, the complications and uncertainties of 
diffusion and mixing are removed from the problem. As a 
result, differences between the experimental data and numeri- 
cal simulations can be attributed either to numerical errors or 
to improperly modeled chemical kinetics. 

An earlier work by Wilson and MacCormack 8 focused on 
the steady ballistic-range experiments and demonstrated the 
physical and numerical modeling requirements for accurate 
computations of these flows. It was shown that a primary 
requirement is adequate resolution in the induction zone be- 
tween the bow shock and the energy-release front. Fine grid 
resolution is needed because the mass fractions of the impor- 
tant radical species change exponentially in this region. Adap- 
tive grid techniques offer one way to efficiently distribute 
points; however, such techniques are not always adequate. The 
induction zone may cover such a large region of the flowfield 
that grid requirements still remain large. Furthermore, the 
small grid spacings in the adapted regions may result in unde- 
sired time step limitations. Additional difficulty is added by 
the presence of the complex flow structures seen in the un- 
steady ballistic-range cases because there is no primary direc- 
tion for adapting. 

As an alternative to adapting the grid, the simulations pre- 
sented in this paper use the logarithmic transformation of the 


species conservation equations presented by Sussman and 
Wilson. 9 Other researchers have realized the potential advan- 
tages of such an approach because the exponential growth of 
the radical mass fractions that occurs in the induction zone can 
be accurately represented with fewer points if a logarithmic 
transformation is used. The formulation of Sussman and 
Wilson is unique because it can be implemented in modern 
shock-capturing computational fluid dynamics (CFD) schemes. 

The use of the logarithmic transformation allows accurate 
simulations of both steady and unsteady ballistic-range exper- 
iments on grids more coarse than have been used in prior 
works. 8 * 10 The reduction in grid size, and thus computational 
resources, permits the use of a detailed chemical model, which 
enables direct comparisons of the simulations to the experi- 
ments to be made for the first time. In particular, measured 
oscillation frequencies are compared with those predicted by 
the computations. 

The simulations in this work add insight into the combustion 
instabilities observed in ballistic-range experiments and allow 
previously proposed mechanisms for the instabilities to be as- 
sessed. In particular, this work considers the wave interaction 
mechanism of McVey and Toong 6 and the modified mecha- 
nism proposed by Matsuo and Fujiwara. 10 McVey and Toong 
developed their mechanism by using experimental data and 
analytical calculations. Matsuo and Fujiwara used modern 
computational fluid dynamics techniques, fine grids, and a 
two-equation, two-step combustion model to do qualitative 
simulations that led them to propose their modified mecha- 
nism. In addition to adding greater physical understanding, 
the current method also provides a tool to validate and possi- 
bly “tune” proposed chemical kinetic models. 

The chemical kinetic model used in this paper can be found 
in Wilson and MacCormack. 8 It is the mechanism proposed by 
Jachimowski 11 containing 13 chemical species and 33 reactions 
except that one of the reaction rate expressions (the expression 
most influential for ignition times) has been replaced by an 
expression recommended by another source. This modifica- 
tion was found by Wilson and MacCormack 8 to give better 
agreement with experiment. Additional details on the chemical 
kinetics are found in Wilson. 12 

Logarithmic Form of the Species Conservation Equation 

This section gives an overview of the logarithmic transfor- 
mation of the species conservation equations as presented by 
Sussman and Wilson. 9 The transformation is based on the 
definition — 





= pHCs) 


( 1 ) 


where p is the total density, p s is the density of species s, and 
c s is the mass fraction of species s. A special property of the 
relationship in Eq. (1) is that it can be used along with the total 
mass conservation equation, 


d „ 

- 1 = 17 = 7 ^ 1 . - J > = 


( 2 ) 


to transform the. species mass conservation equation (neglect- 
ing diffusion), 


a d 

dt Ps + dXj 


PsUj = H'j 


(3) 


into a transformed species equation, 


£ 

bt 





P 

TT S Uj = Ws — 
Ps 


(4) 


where Uj is the j th component of velocity, and w, is the chem- 
ical source term for species s . The transformation from Eq. (2) 
to Eq. (4) retains conservation form. 
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In theory, all of the species mass conservation equations can 
be replaced with a corresponding transformed species equation 
in the form of Eq. (4). In practice, however, this is not done. 
Although Eq. (4) is in conservation form, its use does not 
guarantee conservation of mass or conservation of the individ- 
ual elements. Enforcement of these constraints is achieved by 
including element mass conservation equations in the form 

d * 5 * 

( 5 , 

where p* is the total density of element e contained in all of the 
species (the * is used to avoid confusion between element and 
species densities; pf5 ^ p N ). There are as many of these equa- 
tions as there are elements in the chemical model. The total 
density is given by 

N dements 

p= E P* (6) 

e= I 

The chemical reaction mechanism used in this work involves 
13 species: N 2 , 0 2 , H 2 , NO, OH, N0 2 , HNO, H0 2l H 2 0, 
H 2 0 2 , N, O, and H. Since three elements are present in the 
mechanism (N, O, and H), three equations of the form of 
Eq. (5) are required. The densities pj, pj, and p£ are defined 
using expressions that are the sum of the contributions from 
each of the species containing the element of interest. For 
example, the expression for pj$ depends on the densities of N 2 , 
NO, N0 2 , HNO, and N and is written 

Pn - 2 - N * 4- ^ N0 i | pHN0 + Pn 
A/n A/n 2 A/no A/ N q 2 A/ HN o A/ n 

where M s is the molecular mass of species s. The choice of 
which three species conservation equations to replace with 
element mass conservation equations is not unique. The selec- 
tion of N 2 , 0 2 , and H 2 seems appropriate because these species 
do not exhibit exponential mass fraction growth. The remain- 
ing 10 species conservation equations (for species NO, OH, 
N0 2 , HNO, H0 2 , H 2 0, H 2 0 2j N, O, and H) are written in the 
form of Eq. (4). The densities of N 2 , 0 2 , and H 2 are computed 
from relations of the type of Eq. (7). 

Numerical Formulation 

This section presents the Euler equations with the logarith- 
mic form of the species conservation equations. The 13 species 


• Experimental Shock Position 
O Experimental Energy Release Front 




Fig. 3 Density contours from a numerical simulation of the Mach 
6.46 experiment in Fig. 1: a) flux-vector splitting, b) upwind TVD. 


mass conservation equations of the conventional formulation 
have been replaced by equations of similar form except that 3 
are elemental mass conservation equations [see Eq. (5)] and 10 
are in the logarithmic variable tt s [see Eq. (4)]. The momentum 
and energy equations in the logarithmic formulation are the 
same as they were in the conventional formulation. In vector 
notation the governing equations are written 


dU dF dG 
dt dx dy 


= W 


( 8 ) 


where U is the state vector, F and G are the convective (invis- 
cid) flux vectors in the x and y coordinate directions, respec- 
tively, and IF is the vector of source terms. These vectors are 
given next: 
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where E v is vibrational energy per unit volume, w v is the source 
term representing vibrational relaxation, E is the total energy 
per unit volume, and p is pressure. A separate vibrational 
energy equation is not required for the present computations 
but has been retained from previous work. For the cases 
presented here, thermal equilibrium is simulated by using a 
small time constant for vibrational relaxation. _ 

In Eq. (8), pressure is a homogeneous" function of degree 
one in the conserved variable s T ^and the flux Jacobians (dF/ 
bU and dG/dU) have the s arr^ ^genval ues as that of the 
formulation with the convenrign^forin of the species mass 
conservation eq uation ff ^Therefo re, existing numerical tech- 
niques can be applied, in a 'st raightforward manner. In this 
work, two comrhohlyfusSl shock-capturing numerical tech- 
niques are empl(^K^els a modified Steger-Warming flux- 
vector splitting technique based on work by MacCormack 13 
and Candler 14 (this method shall be referred to as the flux-vec- 
tor splitting scheme) and the other is the Harten-Yee upwind 
total variation diminishing (TVD) (non-MUSCL) scheme 15 
(referred to as the upwind TVD scheme). Both techniques are 
spatially second order and solved in a point-implicit manner. 
That is, the convective terms are treated explicitly whereas the 
source terms are treated implicitly in time (because of the 
small time scales of the combustion chemistry and vibrational 
relaxation). Therefore, a block inversion is required at each 
point at each time step. Details about the flux -vector splitting 
scheme using the logarithmic formulation can be found in 
Wilson. 12 
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Several practical details of the logarithmic formulation are 
addressed next. Because zero mass fractions are not allowed in 
the transformed equations, the mass fractions of the radical 
species must be set to small arbitrary values in the freestream 
(10" 10 for this work). One-dimensional calculations show that 
the solutions are relatively insensitive to the freestream radical 
mass fractions, as long as they remain sufficiently small. An 
advantage of using logarithms is that mass fractions remain 
positive. Limiters are required to prevent the mass fractions of 
species not represented by a logarithm from becoming negative 
(N 2 , 0 2 , and H 2 ). In practice, the logarithmic formulation is no 
less robust than the conventional formulation for the flows 
presented here. The cost per time step is increased by 10-20% 
over the conventional formulation due to the cost of comput- 
ing logarithms and exponentials and also due to the increased 
complexity of the Jacobian of the chemical source terms. 

The logarithmic formulation has limitations at interfaces 
(contact discontinuities) separating gases of different chemical 
composition. The inherent numerical diffusion for each spe- 
cies across such an interface is modified by the logarithmic 
transformation in such a way that the interface is distorted. 
The use of the transformation is advantageous in this work 
because the effect of this distortion is small. In simulations of 
other flows (e.g., the interface in a shock tube), such distor- 
tions may not be tolerable, and the use of the logarithmic 
transformation may not be warranted. 

Numerical Simulations 

This section presents numerical simulations of ballistic- 
range shadowgraphs from experiments by Lehr. 4 * 5 The partic- 
ular experiments of interest used spherical-nosed projectiles 
with cylindrical afterbodies of 15-mm diameter. The cases in- 
clude a range of Mach numbers so that both steady and un- 
steady flows are represented. All of the cases to be considered 
have freestream with a temperature of 292 K, a pressure of 320 
mm Hg, and a premixed stoichiometric hydrogen-air mixture. 
At these conditions, the detonation wave speed was reported 
by Lehr to be Mach 5.11. 

Before the unsteady cases are considered, a simulation of the 
steady Mach 6.46 case in Fig. 1 is shown to demonstrate the 
advantage of the logarithmic formulation and to give an un- 
derstanding of why it is practical to use this formulation for 
the unsteady cases. Figure 3 contains density contours for 
computations of the Mach 6.46 case using flux-vector splitting 
and the upwind TVD methods, respectively. Both computa- 
tions use the 52x52 grid shown in Fig. 4. The experimental 
bow shock and energy-release front positions are overplotted 
on the density contours and show that both numerical methods 
provide good agreement with the experiment. A numerical 
simulation using the conventional form of the species conser- 



Fig. 4 52x52 grid used for the calculations presented in Fig. 3. 



Fig. 5 Density contours from a numerical simulation of the Mach 
4.79 experiment in Fig. 2: a) flux-vector splitting, b) upwind TVD. 



Fig. 6 Computed shadowgraphs from a numerical simulation of the 
Mach 4.79 experiment in Fig. 2: a) flux-vector splitting, b) upwind 
TVD (shadowgraphs computed by Yates 16 ). 


vation equations can be found in Wilson and MacCormack. 8 
Good agreement with experiment was found in that work as 
well, but it required an adapted 321 x64 mesh (321 points 
around the body and 64 points normal to the body). The grid 
for the calculation using the logarithmic formulation is more 
than eight times smaller because the induction zone is ade- 
quately resolved with fewer grid - 

Simulations of Lehr’s unsteady Mach 4.79 case shown in 
Fig. 2 are now presented.^Density contour plots for both the 
flux-vector splitting andjupwind TVD methods are presented 
in Figs. 5a and 5b, ? respectively. The figures represent one 
point in time and ope point in the period of the instability (the 
two solutions afe riot at exactly corresponding points within a 
period) and show that both computational methods predict 
unsteady behavior. As in the steady computation of Fig. 3, the 
density contours exhibit an outer bow shock followed by an 
induction zone and an energy-release front. The most obvious 
difference between the unsteady calculation and the steady one 
is the pulsing structure of the energy-release front. These 
pulses are seen in the experimental shadowgraph with the off- 
axis view in Fig. 2b. The computations were started from a 
fiowfield that was initially set to freestream conditions every- 
where and advanced in time without chemistry until a bow 
shock was established, at which time the combustion chemistry 
was turned on. No other special procedures were used. The 
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Fig. 7 Enlargement of the Mach 4.79 ballistic-range shadowgraph in 
Fig. 2a (courtesy of H. F. Lehr). 


Compression Waves 



Fig. 8 Schematic diagram of the periodic flowfield structures seen in 
shadowgraphs of ballistic-range experiments. 


wise separation between the vertical lines in the computed 
shadowgraphs compared with the experiment, and thus the 
computations predict a lower frequency than the experiment. 
Both numerical methods predict an instability frequency of 
approximately 530 kHz compared with the measured value of 
720 kHz. The primary difference between the flux-vector split- 
ting and upwind TVD solutions is the resolution of the flow 
structures near the shoulder region of the projectile. The flux- 
vector splitting smears these features to a larger degree. 

Further support of the computations is provided by examin- 
ing the experimental and computed shadowgraphs in more 
detail. A schematic diagram of the flowfield that labels the 
different structures is found in Fig. 8. An analysis of these 
structures by McVey and Toong 6 identified the waves extend- 
ing from each pulse of the energy-release front to the bow 
shock as contact discontinuities (McVey and Toong refer to 
them as entropy waves). Their identity was established by ob- 
serving that they are convected with the local fluid speed. 
Compression or expansion waves would move faster than the 
local flow with a relative speed that depends on the speed of 
sound in the mixture. The waves that look like reflections of 
the contact discontinuities off the bow shock do move at a 
speed greater than the local flow and were identified as com- 
pression waves. The dark-light shadings across the waves sup- 
port their asserted identities. All of these features can be ob- 
served in both the' experimental and simulated _shadowgraphs L 

Investigation of other flowfield quantities adds information 
that cannot be deduced from the shadowgraphs. Pressure con- 
tours reveal that each pulse of the energy- release front creates 
a compression wave that travels both upstream toward the bow 
shock and downstream toward the projectile body. Figure 9 
shows the pressure contours at one point in time. The wave 
patterns between the bow shock and the projectile appear quite 
complicated but can be understood by noting that there are 
primarily two families. One family is made up of compression 
waves originating directly from the energy-release front and 
the other family is made up of compression waves from the 
first family that have reflected off the projectile body. 

Instability Mechanism 

The mechanism that causes and sustains the combustion 
instabilities observed in the ballistic-range has not been ob- 
served experimentally. Therefore, plausible explanations for 
the unsteadiness have been developed by extrapolating data 
from outside the nose region (mostly from the shadowgraphs) 
and by transferring knowledge from other flow situations. It 
was not until the work of Toong and his associates that plau- 
sible detailed mechanisms for the regular 6 and large-distur- 
bance 7 regimes were proposed. 


computations used a nearly equally spaced 375x 161 grid. 
Computations on a 375x81 grid yielded similar overall flow 
features and an oscillation frequency near that predicted by the 
375x 161 grid. 12 (Note that the grid required to resolve the 
complex unsteady flow structures negates some of the advan- 
tage that the logarithmic transformation provides. The useful- 
ness of the transformation is the relative grid independence of 
the induction zone length and the instability frequency. This 
allows practical grid refinement studies to be done and adds 
confidence in the accuracy of the predicted frequencies.) The 
code uses approximately 2.8 x 1CT 4 CPU s per time step per 
grid point on a Cray Y-MP. Depending on the grid, the simu- 
lations required between 10,000 and 20,000 iterations and be- 
tween 16 and 70 h of CPU time. 

Figures 6a and 6b contain shadowgraphs computed from the 
flux- vector splitting and upwind TVD flowfields, respec- 
tively. 16 A comparison of these computed shadowgraphs with 
the enlargement of Fig. 2a presented in Fig. 7 clearly shows 
many similarities. Note that the pulses in the energy-release 
front create the vertical line pattern when the axisymmetric 
flowfield is projected onto a plane. There is a larger stream- 
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■, .4 The work of McVey and Toong 6 is discussed further here 
t because it is supported by the present numerical simulations. 
McVey and Toong used what they called a wave-interaction 
model to explain the oscillations observed in the regular 
regime. All of the important processes in this model occur 
between the energy-release front and the bow shock in the 
subsonic region ahead of the projectile nose. As an aid to 
understanding the model, several of the primary components 
of the model are isolated and described first. 

A major component in the McVey and Toong model is the 
interaction between the bow shock and a compression wave. 
Figure 10 depicts a one-dimensional flow at several different 
times (time proceeds from bottom to top). Initially there is a 
stationary normal shock and a compression wave downstream 
of the shock moving toward the shock. At a later time the 
compression wave overtakes the shock, strengthening it and 
causing the shock to move forward. The strengthened shock 
causes a change in the fluid properties behind it. Most impor- 
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Fig. 12 McVey and Toong 6 wave interaction model. 
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Fig. 11 Schematic diagram of a bow shock/compression wave inter- 
action with shock-induced combustion; a) series of one-dimensional 
diagrams, b) x-t diagram. 


tant for the instability is a higher fluid temperature. This 
higher temperature fluid is separated from the fluid that has 
crossed the shock at the original strength by a contact discon- 
tinuity that convects downstream. Additionally, a rarefaction 
wave is created that travels downstream. 

Interesting flow features are created when the interaction of 
Fig. 10 is combined with the flow of a combustible mixture. 
This is the situation presented in Fig. 11a. As in Fig. 10, there 
is initially a stationary normal shock wave with a compression 
moving toward it from the downstream side. Now there is also 
burned fluid downstream of the shock at some induction 
length. As before, the compression wave overtakes the bow 
shock, strengthening it and creating a contact discontinuity 
and a rarefaction wave. Because the fluid on the upstream side 
of the contact discontinuity is hotter, it has a shorter induction 
time and burns more quickly. The result is that for a time there 
are two regions of burning, one at an induction length corre- 
sponding to the original postshock conditions, and the other at 
the shorter induction length corresponding to the strengthened 
shock. Figure 1 lb is anx-f diagram corresponding to the inter- 
action in Fig. 11a. The general features of this diagram are 
recognizable in the McVey and Toong mechanism presented 
next. 

The McVey and Toong mechanism is schematically depicted 
in the x-t diagram in Fig. 12. The diagram shows time evolu- 
tion of the features along the stagnation streamline and is 
explained in the following four steps. ~ 

1) A compression wave on the downstream side of the bow 
shock overtakes the bow shock, causing the bow shock to 
move forward, thus creating a reflected rarefaction (which is 
weak and is ignored in thejnodM|and a contact discontinuity. 
The gas on the upstream fsid eof the contact discontinuity is 
hotter than on the down stream . side because it has passed 
through the strengthened Bow shock. 

2) The rarefaction ? wave propagating upstream from the 
energy-release front overtakes the bow shock, weakening it 
and restoring it to its original strength. The origin of the rar- 
efaction wave is discussed in step 4. 

3) The contact discontinuity created in step 1 convects 
downstream. At a point between the bow shock and the en- 
ergy-release front, a new energy-release front is produced by 
the hotter gas on the upstream side of the contact discontinu- 
ity. The creation of a new zone of combustion, in turn, creates 
compression waves that propagate both upstream toward the 
bow shock (this is the wave that overtakes the bow shock in 
step 1) and downstream toward the projectile. It can be shown 
that the upstream and downstream compression waves are 
necessary to match the fluid velocity jumps across the new 
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Fig. 13 Computed x-t diagrams of density and pressure along the 
stagnation streamline for the Mach 4.79 case using flux-vector split- 
ting: a) density contours, b) pressure contours. 


energy-release front and satisfy conservation of momentum 
(see Alpert and Toong 7 ). 

4) The contact discontinuity (with burning on the upstream 
side) eventually reaches the position of the original energy-re- 
lease front. This extinguishes the original energy- release front 
and creates rarefaction waves that propagate in the upstream 
and downstream directions. The energy-release front then be- 
gins to recede toward the location of the original front as 
colder gas (through the bow shock of original strength) reaches 
the front. 

From Fig. 12 it is seen that the period of an oscillation is the 
time required for the contact discontinuity to convect from the 
shock to the reaction front and then for an acoustic wave to 
travel back to the bow shock. If the postshock conditions are 
known, an analytical expression for the period of observed 
experimental oscillations can be found. Given a period of os- 
cillation provided from an experiment, McVey and Toong 
were able to use the expression to predict the induction time 
for the gas mixture. Their predicted induction times agreed 
with previously published data. Figure 13 contains computed 
x-t diagrams of density and pressure along the stagnation 
streamline for the Mach 4.79 case (the flux-vector splitting 
algorithm was used to produce all of the computed x-t dia- 
grams in this paper), A comparison of these diagrams with the 
McVey and Toong mechanism in Fig. 12 reveals many similar- 
ities. The density contours show a changing position of energy- 
release front in time that is identical to the pattern predicted by 
the wave interaction model (note the jagged edge of the energy- 
release front that appears in both Figs. 12 and 13a). In the 
region between the bow shock and the energy-release front, 
density contours also show that a contact discontinuity is cre- 
ated at the bow shock when a compression wave originating 
from the energy-release front overtakes it. The hotter gas on 
the upstream side of the contact discontinuity reduces the igni- 
tion delay time and causes a new energy-release front to be 
created. The new energy-release front is the source of the 
compression waves traveling in both the upstream and down- 
stream directions. These compression waves are clearly seen in 
the pressure contours of Fig. 13b. The identity of the contact 
discontinuities labeled in Fig. 13a is verified by their absence in 
the pressure contours of Fig. 13b because a contact discontinu- 
ity has no pressure jump across it. 

A phenomenon seen in the computations that was appar- 
ently not anticipated by McVey and Toong is the existence 
and/ or importance of the compression waves reflecting off the 
projectile nose. Figure 13b shows that a compression wave 
created at the energy-release front travels toward the projectile 
nose, reflects, and eventually overtakes the bow shock. In this 
Particular calculation, the compression wave reflects off the 
nose and moves toward the bow shock, overtaking it at nearly 
1 e same time that the most recently created compression wave 
arrives at the bow shock. The coordination between the com- 


pression wave and the reflected compression wave is not al- 
ways exact, and therefore the two waves sometimes hit the bow 
shock at slightly different times. This tends to deform some of 
the structures. The importance of the reflection of the com- 
pression waves off the projectile has been recognized in the 
recent work of Matsuo and Fujiwara. 10 

Figure 14 contains computed x-t diagrams of density and 
pressure from a case at Mach 5.04. As in the Mach 4.79 case, 
the interactions between the bow shock and energy-release 
front correspond to the model proposed by McVey and Toong. 
This case is interesting because the interaction of the wave 
reflected off the projectile is significantly different than that 
seen in the Mach 4.79 case. During the time that a compression 
wave travels from the energy release front to the projectile 
body and back, two new compression waves are created (i.e., 
two periods have passed). The computed frequency is 820 kHz 
whereas the measured frequency is 1.04 MHz. A thorough 
discussion of this case is given in Wilson 12 along with an addi- 
tional simulation of an experiment at Mach 4.18 with a much 
lower oscillation frequency. 

Frequency Sensitivity to the Hydrogen-Air 
Reaction Mechanism 

The nearly constant values of the instability frequency pre- 
dicted by the numerical simulations using different grid sizes 
(375 x 321, 375 x 161, and 375 x 81) and two different numeri- 
cal schemes suggest that the underprediction of the frequency 
in the computations is not caused by numerical error but by the 
chemical kinetic model. Wilson and MacCormack 8 established 
that the blunt body exothermic flowfields are quite sensitive to 
the following chain-branching reaction: 


h + o 2 -oh + o 


In that work, flowfields were computed with two different rate 
expressions for this reaction. One was the original expression 
from the Jachimowski 11 reaction model, and the other was 
from Warnatz. 17 The rate expression recommended by War- 
natz was adopted because it gave better agreement with exper- 
iments for the steady cases (this is the rate expression that has 
been used for all cases presented in this paper). To investigate 
the influence of the chemical kinetic mechanism on the fre- 
quency of the oscillations, the Mach 4.79 simulation is re- 
peated using the Jachimowski rate expression. This rate ex- 
pression gives shorter ignition times at high temperatures 
{T> 1 400 K) than the rate constant recommended by Wamatz. 
This change should lead to an increase in the oscillation fre- 
quency because the induction zone is smaller and therefore the 
travel time between the energy-release front and the bow shock 
is reduced. Consistent with this expectation, the computed 
frequency of the oscillations increased from approximately 
530 to 820 kHz. Since the experimentally determined fre- 
quency is 720 kHz, it is concluded that the uncertainties in the 



Fig. 14 Computed x-t diagrams of density and pressure along the 
stagnation streamline for the Maeh 5.04 case using flux-vector split- 
ting: a) density contours, b) pressure contours. 
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rate constants for the reaction mechanism could explain the 
differences between the experiments and the computations. 
The sensitivity of the oscillation frequency to the chemical 
reactions make the simulations a useful validation tool for 
kinetic models. The unsteady cases appear to provide a better 
measure of a kinetic model than the steady cases because the 
oscillation frequency can be determined with greater precision 
than the positions of the bow shock and energy-release front. 

Finally, note that the arrival times of new compression 
waves and reflected pressure waves at the bow shock some- 
times differ significantly (this occurred in the Mach 4.79 case 
with the Jachimowski rate expression 12 ). The separate interac- 
tions between the bow shock and each wave lead to separate 
contact discontinuities and cause the x-t diagrams to appear 
less organized. Further investigations are required to assess the 
precise contribution to the combustion instability of the com- 
pression waves reflected off the projectile nose. In the present 
calculations they appear to have a secondary effect, but it is 
possible that they are more important to the instabilities seen 
in the large-disturbance regime investigated by Alpert and 
Toong. 7 


Conclusion 

The successful simulations of both steady and unsteady 
exothermic ballistic-range experiments have added to the un- 
derstanding of shock-induced combustion and have led to im- 
provements in the numerical techniques for these types of 
flows; namely, proper grid resolution is a primary challenge of 
simulations with combustion phenomena, and these grid re- 
quirements are reduced by the use of the logarithmic transfor- 
mation of the species conservation equations. The accuracy of 
the simulations is supported by grid refinement studies and the 
use of two different numerical methods, a flux-vector splitting 
scheme and an upwind TVD scheme. Further validation is 
provided by comparison to past analytical and numerical pre- 
dictions of the instabilities. Finally, the simulations have reaf- 
firmed the value of ballistic-range experiments as a source of 
data for high-speed flows with combustion and have demon- 
strated their usefulness for validating numerical methods and 
chemical kinetic mechanisms. It is shown that differences in 
the experimental and computed instability frequencies are 
probably due to uncertainties in the chemical reaction rates 
and therefore the experiments provide a means of validating 
proposed chemical kinetic models. 

Acknowledgments 

Funds for the support of this study have been allocated by 
the NASA Langley Research Center and by the NASA Ames 
Research Center under the joint research interchange number 
NCA2-455 and are greatly appreciated. This material is also 
based on work supported by a National Science Foundation 
Graduate Fellowship and by NASA under Hypersonic Train- 
ing and Research Grant NAGW 965. The authors would also 
like to acknowledge the computer resources provided by the 
Aerothermodynamics Branch of NASA Ames, the Theoretical 


Flow Physics Branch of NASA Langley, and the National 
Aerodynamic Simulation (NAS) Program. We are most grate- 
ful to H. F. Lehr of Institute Saint Louis for providing original 
photographs of his ballistic-range experiments and to Leslie A. 
Yates for computing shadowgraph images from our flowfield 
solutions. We also acknowledge the support and encourage- 
ment of our advisor at Stanford University, Robert W. Mac- 
Cormack. 


References 


JRuegg, F. W., and Dorsey, W. W., “A Missile Technique for the 
Study of Detonation Waves,” Journal of Research of the National 
Bureau of Standards , Vol. 66C, No. I, 1962, pp. 51-58. 

2 Behrens, H., Struth, W., and Wecken, F., “Studies of Hyperveloc- 
ity Firings into Mixtures of Hydrogen with Air or with Oxygen,” 
Tenth Symposium (International) on Combustion, Combustion Insti- 
tute, Pittsburgh, PA, 1965, pp. 245-252. 

3 Chernyi, G. G., “Onset of Oscillation in the Presence of Detona- 
tion Wave Weakening,” Journal of Applied Mathematics and Me- 
chanics ; Vol. 33, No. 3, 1969, pp. 451-461. 

4 Lehr, H. F., Institute Saint Louis (ISL), Final Rept. 20/71, Saint- 
Louis, France, 1971. 

5 Lehr, H. F., “Experiments on Shock-Induced Combustion,” As - 
tronautica Acta, Vol. 17, Nos. 4 & 5, 1972, pp. 589-596. 

6 McVey, B. J., and Toong, T. Y.» “Mechanism of Instabilities of 
Exothermic Hypersonic BIunt-Body Flows,” Combustion Science and 
Technology, Vol. 3, No. 2, 1971, pp. 63-76. 

7 Alpert, R. L., and Toong, T. Y., “Periodicity in Exothermic Hy- 
personic Flows About Blunt Projectiles,” Astronautics Acta, Vol. 17, 
Nos. 4 & 5, 1972, pp. 539-560. 

8 WiIson, G. J., and MacCormack, R. W., “Modeling Supersonic 
Combustion Using a FuIIy-Implicit Numerical Method,” AIAA Jour- 
nal, Vol. 30, No. 4, 1992, pp. 1008-1015; see also AIAA Paper 90- 
2307, July 1990. 

9 Sussman, M. A., and Wilson, G. J., “Computation of Chemically 
Reacting Flow Using a Logarithmic Form of the Species Conservation 
Equations,” Proceedings of the 4th International Symposium on 
Computational Fluid Dynamics (Davis, CA), Vol. 2, Sept. 1991, pp. 
1113-1118. 


10 Matsuo, A., and Fujiwara, T., “Numerical Simulation of Shock- 
Induced Combustion Around an Axisymmetric Blunt Body,” AIAA 
Paper 91-1414, June 1991. 

11 Jachimowski, C. J., "An Analytical Study of the Hydrogen-Air 
Reaction Mechanism with Application to Scramjet Combustion,” 
NASA TP-2791, Feb. 1988. 

12 Wilson, G. J., “Computation of Steady and Unsteady Shock-In- 
duced Combustion over Hypervelocity Blunt Bodies,” Ph.D. Thesis, 
Stanford Univ., Stanford, CA, Dec. 1991. 

I3 MacCormack, R. W., “Current Status of Numerical Solution# of 
the Navier-Stokes Equations,” AIAA Paper 85-0032, Jan. 1985. 

,4 CandIer, G. V., “The Computation of Weakly Ionized Hyper- 
sonic Flows in Thermo-Chemical Nonequilibrium,” Ph.D. Thesis, 
Stanford Univ., Stanford, CA, June 1988. 

I5 Yee, H. C., “A Class of High-Resolution Explicit and Implicit 
Shock-Capturing Methods,” NASA TM 101088, Feb. 1989. 

I6 Yates, L. A., “Images Constructed from Computed Flowfields,” 
AIAA Paper 92-4030, July 1992. 

I7 Warnatz, J., “Rate Coefficients in tfie]C/H/0 System,” Combus- 
tion Chemistry, edited by W. C. Gardiner Jr., Springer-Verlag, New 
York, 1984. Chap. 5. 


it 


OaKHWL PAGE 16 
OF POOR QUALITY 



Hypersonic Aerodynamic Decelerators Design Using CFD 
Part II - Turbulent Computations 
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Abstract 


Aerodynamic and heating predictions for a variety of hypersonic decelerators have been per- 
formed with a number of Computational Fluid Dynamic (Cr D) solvers in order to determine the 
design parameters. Sample results from the computed solutions are presented ior a number of 
design. In addition to thermo-chemical non-equilibrium effects, turbulent effects are significant in 
determining the efficiency of the decelarator design. Efficient approximation of treating the thermo- 
chemical, non- equilibrium gas as an idealized gas has been proven to have significant advantage in 
the decelerates design and the use of well established algebraic turbulence models are explored in 
this study. The accuracy of the computations are evaluated by comparing the flow-field solutions 
and aerodynamic coefficients with experimental measurements. 

Introduction 

The use of supersonic and hypersonic drag-devices or decelerators have been experimentally 
investigated for many years with applications ranging from high-performance re-entry vehicle re- 
covery to booster-pavload collision avoidance systems. Supersonic and hypersonic drag-device con- 
figurations that have been examined in the past include attached and tethered ballute, ' which 
is a combination balloon and parachute, and rigid flare, 7 " 10 among others. Recent advances made 
in modelling hypersonic non-equilibrium flows and flow* solvers have significantly improved *,he 
design of complex drag devices. 11-13 The drag device designed investigated and designed by the 
present authors using the CFD solvers include ballutes and flares 1 " for solid rocket motor (SRM) 
casing earth-capture mission and petal-flares for Lunar (and Mars) return vehicles (LRV) . These 
designs are shown in Figure 1-3. In the case of ballute and flare design for the SRM, the CFD 
solvers were used to determined the size and shape parameters that will produce sufficient drag and 

* 


Research Scientist, Eloret Institute. Mailing Address: NASA Ames Research Center, MS 
230-2, Moffett Field, Ca 94035. AIAA Member. 


1 


survive reentry speeds of 9 Km/sec during reentry. The drag coefficient required were prescribed 
based on the ballistic coefficient for a no-skip reentry trajectory. Computations were performed 
for the LRV decelerator to determine the drag and heating for a variety configurations and shapes. 

At re-entry speeds around 8 km/sec and at freestream conditions corresponding to an altitude 
of 80 km, the flow around the fore-body and the decelerator section remains in thermo-chemical 
non-equilibrium. The flowfield around the decelerator section can be highly complex and often 
a separated flow region is established. In addition to the separated zone, the inviscid- viscous 
interaction may produce either a strong or an oblique shock wave and the corner sIioca interacts 
with the fore-body shock. The hypersonic forebody is not expected to be turbulent, but the flow 
field near corner regions of the decelarator, due to flow separation, can be dominated by turbulence. 
The computed flowfield around the three different decelerators are shown in Figures 4-6. The non- 
equilibrium issues are addressed in Part I , one of two papers proposed to be presented together. 
The effectiveness of the decelerator depends on the strength of the comer shock and the size of the 
separated domain. Turbulent mixing can significantly affect the corner flow domain and hence the 
efficiency of the decelarator. In addition, if local heating is important to the design, validated CFD 
solvers are needed to be used in the design process. The viscous-inviscid interaction is non-linear 
and complex and the effect of the non-equilibrium, turbulent state oi the gas on the inviscid- viscous 
interaction is difficult to quantify. In addition to the flow in the decelarator section, the flow in 
the base region of the decelerator may influence the upstream flow. In the case of the ballute and 
axisymmetric flares, the base flow is not significant. But in the case of the LRV decelerator, the 
base flow interaction with the upstream flow is significant and the complete flow field including 
the base region need to be solved during the design cycle and the turbulent mixing and modelling 
becomes significant. In the case of the LRV decelarator, parr oi the flow is expanded through the 
gap into the base region and base recirculation is set up. Complex vortex patterns are set up due 
to the three-dimensionality oi the flow and the three-dimensional computational predictions are of 
great significance to the design. 

A number of CFD solvers used by the present authors 12 ’ 13 and others 11 >15 46 ’ 1( in comput- 
ing the hvpervelocity decelerators can be classified into two categories, namely thermo-chemical, 
non-equilibrium solvers 1 ^ ,16 and effective real-gas solver. 11 ,1 ~ fl In the thermo-chenncal, non- 
equilibrium solvers, the Navier-Stokes equations including the thermo-chemistry model equations, 
represented by the species and energy conservation equations, are solved. In the effective real-gas 
solver, the gas model is approximated by an idealized gas with an effective specific heat ratio 
approximating the global non-equilibrium state of the gas and the Navier-Stokes equations alone 
are solved. The effective- real- gas solver is computationally more efficient but less accurate com- 
pared to the full thermo-chemical non-equilibrium solver. The three-dimensional non-equilibrium 
solvers are prohibitively expensive requiring 40-60 CPU hours on CRAY-2 for one three-dimensional 
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solution 15 , but the effective-real-gas solvers 12 |17 require very reasonable amount of computer time 
(1-2 CPU hr.). The non-equilibrium issues will be addressed in another paper and in this paper, 
the turbulent model issues, for the idealized-real-gas solver, are explored. 

Numerical Method 


Effective-7 Ideal-Gas Solver 

The use of CFD in the design often requires a number of solution to be computed. These 
solutions are used to evaluate the impact of various design parameters on the aerodynamic and 
heating and to choose the best of the design variables. In order to calculate solutions in a timely 
and Tost effective manner, an ideal-gas solver with an effective specific heat ratio, 7« ( ,i has been 
used in the design process. 12 ,1,J The underlying assumption is that the global effect of a real- gas 
flown eld can be simulated sufficiently accurate , by using an effective value of 7. This technique 
has been used successfully to predict drag and moment coefficients produced by blunt bodies in 
real-gas flows 17 using the effective-7 ideal-gas solver developed by the first author. This solver uses 
an implicit, finite-difference, flux- difference splitting method due to Roe and is suitable for both 
axisymmetric and fully three-dimensional flow problems. The idealized solver is very robust and 
requires only 10 to 15 minutes of CPU time on Cray-2 for an axisymmetric calculation and between 
1 and 2 CPU hours for full 3-D solution. 

To model turbulence, algebraic models 19 were used for wall layer region and a modified mixing 
length model was utilized in the base region. The algebraic models are well suited for simple 
2-dimensional and axisymmetric flows, but require carefull fine tuning for complex 3-dimension<il 
flows. Work is currently in progress to include a compressible 2-equation model/ 0 Full detail 
regarding the the turbulence modelling effort will be given in the unai paper. 


Sample Results and Discussions 

A number of design studies have been performed with the effective- 7 ideal-gas solvers 1 ' by 
the present authors. Some of the computed solutions were shown in Figs. 4-6. The fore-body 
region in these cases is dominated by the non-equilibrium effects and the flow remains laminar. 
The computational studies performed with the effective-7 ideal-gas solver 1 ' were validated by com- 
paring the computed results with ballistic range measurements of drag, lift and moment coefficient 
for a three-dimensional Aero- Orbital Transfer Vehicle geometry. These comparisons prove that 
the effective-7 approach was reasonably accurate in predicting aerodynamic forces at moderate 
hypersonic speeds for blunt fore-body flows. 

The aerodynamic decelerators, in addition to the blunt forebody, have the added complexity 
of the corner flows, where turbulent mixing will affect the flow field. To validate the present solver, 
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computed solutions to 2-dimensional corner flow problems are compared with experimental results. 
A laminar 2-d corner flow comparison is shown in Fig. 7. In this experiment, the flow upstream 
of the corner remained laminar and through the separation bubble, the flow transitioned from 
laminar to turbulent. The heat transfer comparison shown in Fig. 7 agrees well in the laminar 
region and as expected, the difference between the computations and the experiments are in the 
turbulent region. Computed results using a a number of turbulent models and comparisons will 
be presented in the final paper. 

In addiiton to validating the solver for the corner region, the solution accuracy is need to be 
assesed for flows with complex 3-d effects. The computed LRV deceiarator ftownled solution are 
shown in Fig. 8. The flgure on the top left corner shows the simulated surface oil-flow patterns 
along" with the surface pressure. The vortical, recirculating flow near the compression-corner region 
can be seen to be highly complex. The invscid-viscous interaction effects can be visualized from 
the figure on the top right corner. In this figure, the particle traces along the symmetry plane 
and the shape of the outer shock surface are shown. The compression corner region produce 
corner shocks and this shock interact with the outer shock. Tne corner shock, in addition to 
raising the surface pressure, causes local heating and the local heating levels are of paramount 
significance to the design. Only a limited amount of experimental results are available for similar 
flows. Computations are underway to simulate and compare with these results. The experimental 
drag coefficient 14 is shown in Figure 9. The effect of the gap ratio on the drag coefficient is seen 
to be non-linear and 3-d and turbulent effects are beleived to be the cause. Validation efforts and 
final comparisonswill be reported in the final paper. In addition, design efforts earned out for the 
LRV deceiarator will be discussed. 


Summary 


The use of CFD in the design of hypervelocitv accelerators have proven to be very attractive. 
In the final paper, turbulent computations for a number of decslerator design will be reported 
and the preliminary results indicate that the effective-*/ approach provides a reasonably accurate 
computational procedure for the design of complex, hyperveiocity decelerators. 
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Figure 6. Pressure Contours - Petal Decelerator (symmetry plane) 




Figure 7. Stanton Number Comparison for 2-D Corner Flow Problem 
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Figure 9. The Effect of Gap Ratio on the Drag Coefficient (Ref. ) 
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Abstract 


Comparison is made between the results obtained from a state-of-the-art flow and 
radiative model and bow shock vacuum ultraviolet (VUV) data obtained the t recent Bow 
Shock 2 night Experiment. An extensive data set was obtatned from onboard rocket 
measurements at a reentry speed of 5km/sec between the altitudes of approx, irately 65- 
85km. The data, a description of the NO photoionizarion cell used, and interpretation o 
the data will be presented. The primary purpose of the analyses is to access the utility o 
the data in the context of a radiation model appropriate (o the flight conditions o ow 
J Shock 2. Theoretical predictions based on flow modelling discussed earlier and a new 

radiation model are compared with data. 


Introduction 

The first bow shock flight experiment (April 1990) measured die ultraviolet radiation 
from shock-heated gas in the nose region of a 0.1016m nose radius rocket traveling a. 
3 5km/sec at altitudes from 40 to 70km. Papers providing detads on the mstrumentauon 
and analyzes of results* for the ultraviolet sensors have been presented. The second bow 
shock flight experiment (February 1991) provided similar types of ultraviolet data but or 
reentry at 5 Wsec between the altitudes of about 100 to 65km .» In this paper we 
examine vacuum-ultraviolet data (VUV) obtained from the resonance atomtc transition at 
1304A Such data is available only from the second flight. This data has received lower 
priority titan that of tire ultraviolet instruments; however, as will be shown herdds cructa! to 
understanding the shocklayer modelling in its entirety. The reduction of this data was 
complicated by rocket precession during reentry caused by incomplete second g 
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separation. Use of onboard navigational data has permitted us to correlate sensor - 
stagnation streamline angle as a function of time. 

For early time portions of the reentry the atomic oxygen concentration in the shocklayer 
may be sufficiently low such that the resonance radiation is optically thin. In this case the 
radiation obtained is directly proportional to the amount of atomic oxygen tn the flow The 
degree to which theory can match the density dependence of the data indicates the fidelity o 
die fiow thermochemical modelling. For example, a.omic oxygen colliding with N 2 is the 
precursor process to the formation of NO for the fhgh. conditions of the bow shock : fhg t 
experiments. In Ref. 4 we have discussed die uiTW^eemen, “ theory 
and experiment for the 2300A NO radiation, particularly atalntudes £* 


and experiment lor tne v ' v ^ 

may be due tolMtl&iw^i nhfrediCTrotr 

/#NO concentration, its temperature, oShoice of the appropriate governing temperature 
Tthe radiation model. Of these possibilities the first is considered the most likely source 
of difficulty. Since the kinetic processes of NO and O am closely coupled it is unlikely that 
a correct flow solution will be obtained unless both the 2300 and 1304 radiation data sets 

can be modelled. 

At later time portions during reentry flower altitudes) the atomic oxygen concentration 
may increase so that self absorption becomes significant. For these conditions die 
modelling of the shocklayer spatial O concentration and temperature dependence, especially 
in die boundary layer, becomes important. Unlike the 2300A radiation which depends 
primarily on the spatial portion of the shocklayer corresponding to the peak radiation in e 
flow, these conditions would f urther test the ability to model the boundary lay er^. . — ^ ^ 

Xhe^toti^deT^. has ^ 

Excitation of atomic oxygen to the ’S resonance state is achieved by collies with nHWT~ 
species rather than electrons. There are sufficient numbers of collisions such that a steady 
state assumption for the upper state population is assumed and self absorption of 13MA 
atomic oxygen resonance radiation is treated by classical radiative transport theory. The 
NEQAIR 2 model with certain modifications (which will be enumerated) was the analytic 

tool used to compute resonance radiation. 5 


Experimental Set-O P 

Atomic oxygen 1304A radiation measurements were planned for the first flight 
experiment with an NO-filled, CaFj window ionization chamber located in the stagnation 
region of the rocket. Due to quariz window spectial cut-off, the sensor diode was sealed 
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directly to the inside of the payload dome and viewed the ionized gas through a small hole 
of 0 5mm diameter. The VUV gas diode was compared against an NBS calibrate 
diode to yield an absolute calibration for the 1304A radiation measurement. The total 
dynamic range of this device was limited by unceraindes in die amplifier baseline drift that 
occurred due to the proximity of the electronics with the telemetry transpon er an 
unanticipated shielding requirements. Thus 1304A radiation data obtained from the 
Bow Shock flight is considered unreliable. 

VUV measurements were taken during the second Bow Shock flight with an NO 
photoionization cell in a configuration similar ,0 the firs, flight. Major design changes to 
die preamplifiers from the first flight permitted data to be obtained from 85 to 65 m 
alutudes Figure 1 shows the celt and its geometry relative to the rocket payload axis^ 
Figure 2 shows the spectral response of the CaF 2 window which can be seen to be muc 
broader than the O transitions. Since there is essentially no atomic nitrogen present in the 
shocklayer under these flight conditions the radiation obtained from this measurement is 

assumed to be solely that due to atomic oxygen. 
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Fig. 1. NO photoionization cellusetj to measure atomic oxygen radiation. 
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Data 


Figure 3 shows a summary of the atomic oxygen 1304A radiance data that was 
obtained during the reentry portion of the second flight along with the corresponding 
altitude of the rocket. The atomic oxygen sensor was located 37° from the rocket 
longitudinal axis. The highly oscillatory nature of the signal is due to rocket precesston 
which changes the location of the stagnation point as a function of time dunng the flight 
and is more noticeable in the data collected from off-axis sensors. To permit quantitative 
comparison of data with models, the angle between the sensor and the effective stagnation 
streamline (7) must be known. Using the onboard navigational data supplied by San a 
National Laboratories and transformations of the sensor location into the onboard 
navigational coordinate reference frame* ,7 can be computed as a function of time. Figure 
J 4 shows an example of radiance correlated with 7 during a time region selected at ~7>km. 
The radiance peaks at small 7 (U. closer to the stagnation point) and is reduced at the 
largest values (i.e. towards the wake). The time correlation appears to be good and can 
provide detailed information on the angular variation in the flow. 
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Fig. 3 Atomic Oxygen 


Radiance and Vehicle Altitude as a Function of Time after Launch 
- Reentry Portion of Flight 
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Fig. 4 Comparison of Stagnation Streamline - Sensor angle and VUV Radiance as a 

Function of Time after Launch 
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Ra diation Model 

Modelling of 1304A resonance radiadon has been considered by Whiting and Park 5 to 
evaluate headng rates for a blunt body reentry at lOknr/sec (Aeroassist night Experiment). 
The major excitation processes under AFE conditions are produced by electron collisions 
with O. For the Bow Shock 2 Flight experiment the atomic radiauon observed is due to 
neutral ■ O collisions. If it is assumed that there are enough collisions such that the 
electronic states are in a Boltzmann distribudon at the heavy particle translation temperature 
orders of magnitude more radiation is predicted than was observed. Therefore a rate 
equation approach must be considered. A new model is proposed that accounts for neutral 
excitations with O; but, assumes that the upper state 5 S population can be obtained from the 
steady state solution of the appropriate rate equations. 

Table 1 shows the atomic oxygen energy levels for the first five electronic states. 
Rates for neutral coUisional excitations among the firs, three suites have been measured and 
are high for the temperatures obtained in the shocklayer. Therefore these levels are 
assumed to be in a Boltzmann distribution at the heavy particle translauon temperature The 
strongly allowed dipole resonance 3P-»3 S transition can be modelled as the three level 
system shown in Fig. 5. Collisions of the O('D) state with neutrals will produce a 
transition to the 5 S state, although this rate has probably not been explicitly measured, n 

estimate of that rate will be used based on data given in Ref. 8. 


Table 1. Atomic Oxygen Spectral Properties 


State Number 

1 

2 

3 

4 

5 


Term Symbol 
3p 

IS 
5S 
3S 


Degeneracy 

9 
5 
1 
5 
3 


Energy Level (cm' 
78. 

15,868. 
33,792. 
73,768. 
76,795. 


In the optically thin limit the rate equations for this three level system are 


/ 


= k f I0( 3 S)l[M]-k r [0( 1 D)][M]-k 2{ I0( 1 D)l[M]+k 2r [ 0( 3 P)][M] 


d[Q( S)]_ _k ^0( 3 S)] [M] +k f [0( 1 D)] [M] -A[0( 3 S)] 


( 2 ) 


Under steady state assumptions the derivatives in Eqs. (1) and (2) are set to zero and an 
expression for the upper state ^S concentration is obtained. 
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[0( 3 S)] qss 


[0 foxVQ,nft _ [0?Vf\QJQ?h 


n k + A< 


t f 

[ Q 2. 


(3) 


where n t is the total free stream density of species M, A represents the Einstein coefficient 
of 3.8xl0 8 sec 1 for the resonance transition, Qj represent the partition funcuons of the 
states given in Table 1, and k f is the quenching rate from the 3S to ] D state. It is interesting 
to note that since the Einstein A coefficient is about a factor of 100 greater than that for NO 
2300A radiation, 1304 A radiation is collision limited whereas for NO collisional quenching 


dominates. 


Since the production of atomic oxygen is second order in n t , Eq. (3) predicts that the 
atomic oxygen radiation should be proportional to n t 3. To test this hypothesis the density 
dependence of the data can be examined. Figure 6 shows data taken at altitudes between 
80-70km sorted to include measurements that occurred at an angle of 36±0.2 from the 
stagnation streamline. The oscillations can be further reduced with tighter sorting; 
however, the third law dependence in n t can be seen to give the best fit. 



Fig. 5 Three Level Atomic Oxygen System 


Calculations 

The computational technique that has been used to simulate the stagnation region 
flowfield is that given in Refs. 2 and 4 which discusses enhancements to the original flow 
modelling of Candler and MacCormack.’ Briefly, the shock heated air is modelled by 
eight chemical species (N 2 , 0 2 , NO, NO*. N 2 *. N. O. and e') that arc allowed to react 
with each other at finite rates. Three independent temperatures in die flow are modelled: 
heavy particle translational, molecular rotational, and a combined molecular vibrational and 
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Fig.6 Comparison of Atomic Oxygen Data and Density dependence, 

electron translational temperatures. The atomic oxygen concentration and heavy particle 
temperature obtained at each point in the flow have been used to compute the 1304 A 
radiation. Two flow models have been considered, a thcrmochemical "baseline” version 
and a modified solution based on alternative reaction rales suggested by Wray. 4 These 
solutions differ substantially In the concentration of ground state O predicted between 70- 
80km. At 80km there are factors of 100 difference in atomic O concentration between the 

two solutions, 

Equation (3) has been incorporated into the NEQAIR 2 model. The calculations have 
modelled the dependence of the upper state population on optical thickness through the use 
of the escape coefficient. This aspect of the calculation may be refined further in the final 
version of this paper. The full radiative transport capability of the NEQAIR model has 
/ been used to evaluate 1304 A radiation at altitudes of 85, 80, 75, 71, and 65km with both 
/ flow solutions. A value of k f of 8xl0-iW/sec at 300K obtained from Ref.8 was used. 
Although the actual value of the quenching coefficient still needs to be obtained, n is fairly 
constant over the altitude mnge considered here. Hence if one flow model is supenor to the 
other, the shop* of its predicted radiance should be in better agreement with the data. 

Figure 7 shows the results of the calculations for the full radiative transport radiance at 
the body ( l.e the location of the 1304 sensor). For convenience the corresponding 
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Fig. 7 Comparison of New NEQAIR 2 Calculations with Data - Full Radiative Transport 


altitudes are also indicated. The data points at 85 and 65km should have large error bars 
associated with them. At the highest altitude, the data is becoming shot noise limited 
making it difficult to correlate exactly sensor viewing location as a function of time in the 


precession period. At 65km there remains uncertainties regarding the performance of the 
J photoionization ceU, which still need^to be clarified. These error bars will be specified in 
the final version of this paper. Both flow models predict the shocklayer to be optically thin 
for altitudes^)75km. At 71km the baseline model predicts self absorption in the 
shocklayer to be more important than the Wray flow modelling. Both flow models predict 
the self absorption in the shocklayer to be substantial by 65km. If the comparison between 


theory and data is restricted to altitudes between 70 to 80km where all models predict the 
flow to be optically thin and the data density dependence and radiation modelling are 
consistent Fig. 6) neither calculation pred^ts the correct radiation density 
dependenc^l^eforeftfie Ihermoctfemicaf mocfellrngof atomic ^x^^^^ation 
should be further examined. 


The disagreement with theory and data for the 65km altitude oF orders of magnitude 
may be due to problems in modelling the atomic oxygen concentration and temperature in 
the boundary layer. Since the NO 230nm radiation was optically thin, the radiation 
computed at the body in that case would be less sensitive to errors in the boundary layer 
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modelling. Grid refinement techniques to emphasize boundary layer modelling will be 
discussed in the final version of this paper. 

Conclusion and Additional Ufo rk 

A radiation model has been successfully proposed to describe neutral excitations of 
atomic oxygen that produce 1304 A resonance radiation. The magnitude of the specific 
quenching rate or the ability to calculate it will be discussed further in the final version of 
this paper. Self absorption has been included in the calculations shown here; but, the 
approximate escape coefficient formalism will be replaced with an improved iterative 
treatment of the radiation field in the final version of this paper. 

The difference in density dependence discrepancy that was obtained for the 2nd flight 
comparison between the 0° viewing 230nm photometers and calculations are seen here for 
the atomic oxygen calculations as well. This data is important since a rigorous test of any 
flow solution should include comparison with the 1304 as well as 2300A radiation. Since 
atomic oxygen is a precursor to formation of NO the correct calculation of its concentration 
is crucial in modelling 230nm radiation. 
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Abstract 


Radiation from Nitric Oxide {NO) band systems in the free-stream of NASA Ames 
20 MW arc-jet wind tunnel was measured and computed. The emission measurement was 
conducted with a 0.3 meter McPherson spectrograph using photographic films in the wave- 
length region from 185 nm to 300 nm. The NO band systems consist of the following bands; 
NO3, NOs , NO,, and N0 1 . The analysis of data indicated that only the ground vibrational 
level {v' = 0) of NO - y and NOs band systems are present. 

Excitation temperatures in the free stream of arc-jet flow were deduced by comparing 
the measured spectra with the calculated spectra using the nonequilibrium radiation code 
NEQAIR. The deduced rotational, vibrational, and electronic temperatures are: T r = 6G0 ± 
50° K, T v < 850 ± 50° K, and T e = 7560 ± 340° K, respectively. The results indicated that 
the free stream flow was in thermal nonequilibrium. 

The Three-Temperature Nonequilibrium Nozzle Flow computer code (NOZ3T) was used 
to predict the arc-jet nozzle flow. The code assumes the early portion of the nozzle flow 
to be in equilibrium. The rest was solved assuming three-temperature nonequilibrium flow. 
Comparing the computed temperatures of T r , T v , and T e to those extracted from the expeii- 
mental spectra the following two important parameters were determined; the cross section for 
transfer of vibrational energy to electronic energy, and the ratio of true vibrational relaxation 
rates to Milikan and White vibrational relaxation rates. 


Experiment 

Figure 1 shows the schematic of the experimental set-up. The operating conditions of 
the 20 MW arc heater were the following: chamber pressure 2. 4 a/m, flow rate 0.0864%/sec, 
current 2000.4 and voltage 2600 E. Spectral data were taken in the free-stream flow 14 inches 
downstream of the nozzle exit[l]. 

A McPherson spectrograph model 218 with a camera attachment was used to collect 
the data. The spectrograph has a speed of // 5.3 and 0.3 meter focal length. The set-up 
used is known as Criss-Crossed Czerny-Turner optical system [2]. It is an unconventional 

’Graduate Research Assistant, Stanford University, Stanford, California. Member AIAA. 

2 Research Scientist, Eloret Institute, Palo Alto, California. Member AIAA. 

3 Research Scientist, NASA-Ames Research Center, Moffett Field, California. Associate Fellow AIAA. 
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mounting arrangement where the light beam from the two mirrors and the plane grating are 
crossed in order to obtain the // 5.3 optics. A Joyce Loebel microdensitometer model Mk 111 
CS was used to measure the optical densities recorded on the film. The microdensitometer 
traces were transferred to a computer to be analyzed and compared with the theoretical 
predictions. Detailed description of the experimental set-up, instrumentation, and data 
collection and calibration procedures will be given in the full paper. 


Data Analysis 


The uncalibrated densitometer trace of the second order spectra is shown in Fig. 2. 
The wavelength calibration was performed using commercial mercury (Hg) and argon (Ai) 
spectral lamps. The mercury lamp was used in the 400nm to 650nm spectral range and the 
argon lamp was used in the 650nm to SoOnm spectral range. From the densitometer traces 
of the mercury and argon spectra and from the known locations of their spectral lines, the 
relation between the lamps spectra and the actual free-stream spectra was developed. 

The relative intensity calibration was done in two parts. First the film was calibrated 
for the relative intensity at specific wavelengths. Second the film was calibrated for the 
spectral sensitivity changes due to the film sensitivity, the spectrograph grating, and the 
entire optical set-up. This calibration was done using a natural density step-wedge filter. 
The resulting intensity trace, corrected for the film response is shown in Fig. 3. 

The excitation temperatures (rotational, T r , vibrational, T v , and electronic, T e ) were 
determined by comparing the calibrated experimental and the calculated spectra. ^The cal- 
culated spectrum was obtained using the Nonequilibrium Air Radiation code (NEQAIR) 
[3.4], The upper limit of vibrational temperature was determined by finding the vibrational 
temperature in the calculations that extinguishes the radiation from the v' = 1 level. The 
electronic temperature was deduced from the relative intensities of A 0-, and NOs- The 
temperature determination methodology is iterative in nature and will be discussed in the 

full paper. 


Theory 

The computer code NOZ3T[5,6] was used to predict the free-stream flowfield through the 
convergent-divergent nozzle of the arc-jet wind tunnel. The entrance and the begining section 
of the nozzle were assumed to be at equilibrium. The rest of the nozzle was solved assuming 
three-temperature nonequilibrium flow. The equilibrium thermodynamic conditions required 
at the entrance of the nozzle are: pressure, velocity, estimated temperature , and enthalpy. 
These parameters were computed using Arc Heater Flowfield computer code (ARCFLO) 

[7,8]. _ . . 

With the operating conditions of the arc-jet wind tunnel given in the preceding section, 

the ARCFLO predicted the following center line values at the exit of the arc-column: T - 
12189° I\, V = 226m/ sec, and h = 60 MJ/kg. The corresponding mass averaged values were: 

T = S91S °K, V = 185 m/sec, and h = 22 MJ/kg. 

The results of calculation w r ere compared with the rotational, vibrational, and electionic 
temperatures obtained from the experiment. Tw-o important parameters; the cross-section 
for vibrational-electronic energy exchange between various species, and the ratio of true 
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vibrational relaxation rates to Millikan- White vibrational relaxation rates[9], were adjusted 
to reach to an agreement between the two sets of temperatures. 


Results 

The spectral bands detected by the spectrograph were identified as the NO y and A 0 5 
band systems in second order. The N0 y band system is denoted by A 2 S + -> X IT, where 
A 2 E + is the upper electronic state and X 2 n is the ground state of the Nitric Oxide (A O) 
molecule. The NOs band system is denoted by C 2 n -*■ A' 2 II where C 2 II is the upper 
electronic state and X 2 n is the electronic state of the A O molecule. 

The NOp band system, which is denoted by (J5 2 II — > A 2 II) and lies in the same spectral 
region as the A r 0 7 , was not present in the data. The absence of this band system indicates 
that the B 2 Ti electronic state was not populated to a significant amount, although it is 
located at nearly the same energy level as the A 2 E + state (see Fig. 4). The B fl state 
does have a slightly larger internuclear distance than the A‘S + state, however this will not 
affect their population ratio under equilibrium conditions. This phenomenon suggests that 
the~upper electronic states of the A O molecule was populated in a nonequilibrium manner 

due to a possible selective population mechanism. 

The NO, band system denoted as (£> 2 £ + -» X 2 n) lies very close to the NOs band system 
in the energy potential curves. However the strong bands of the A O t lie in the \ acuum 
Ultra Violet (VUV) region and only few of the weaker bands are in the spectral region 
of our experiment. These weak vibrational bands of A 0, system (vibrational-rotational 
levels (0,5), (0,6), and (0,7)) lie under the much stronger vibrational bands of NO y systems 
(vibrational-rotational levels (0,0), (0,1), and (0,2)) (see Fig. 5). To verify the existence of 
NO, in the free stream flow, further emission measurements are needed in the VUV spectral 
region where NO y band system does not exist. 

The deduced excitation temperatures are T r = 660 ± 50° A, T v < 850 ± 50° A, and 
T e = 7560 ± 340° A”, and the population ratio of the C 2 n state to the A 2 E state of the 
Nitric Oxide (NO) molecule is 0.42 ± 0.04. The study showes that the cross-section for 
vibrational-electronic energy exchange is in the order of 10 _18 c7n, where as the dividing factor 
to Millikan- White relaxation factor is at least in the order of 10 s . The enthalpy deduced 
from this work is ZSMJ/kg. Figure 6 is a plot of the computed excitation temperatures 

using NOZ3T. 


3 



References (partial list for the abstract) 

[ll Gopaul N.K.J.M., “Spectral Measurement of Nonequilibrium Nitric Oxide in the Free- 
Stream of an Arc-Jet Flow,” Degree of Engineer Thesis, Stanford University, Stanfored, 

1992. 

[2] Mcpherson Instruments, “Instruction Manual for Model 218, 0.3 Meter Combination 
Scanning Monochromator and Spectrograph,” Acton, MA. USA. 

[3] Park, C-, “Nonequilibrium Air Radiation (NEQAIR) Program: User’s Manual, NASA 
Technical Memorandum 86707, 1985. 

[4] Moreau, S., Laux, C., Chapman, D.R., MacCormack, R.W., “A More Accurate Nonequi- 
librium Air Radiation Code: NEQAIR Second Generation,” AIAA Paper 92-2968, 1992. 

[5] Park, C. and Babikian, D.S., “Users Manual for NOZ3T The Three-Temperature Nonequi- 
librium Nozzle Flow Code,” 1992. __ m 

[6] Seung-Ho, L. and Park, C., “Validation of Three- Temperature Nozzle Flow Code NOZ3T,” 
Proposed Paper for the AIAA 28th Thermophysics Conference, Orlando, Florida, 1993. 

[7] Watson, V.R. and Pegot, E.B., “Numerical Calculations for the Characteristics of Gas 
Flowing Axially Through a Constricted Arc,” NASA TN D-4042, 1967. 

[8] Nicolet, W.E., Shepard, C.E., Clark, ICC., Balakrishnan, A., Kesselring, J.P., Suchs- 
land, K.E., and Reese, J.J., “Methods for the Analysis of High-Pressure, high Enthalp} Arc 
Heaters,” AIAA paper 75-704, AIAA 10th Thermophysics Conference, Denver, Colorado, 

1.975. 

[9] Millikan, R.C., and White, D.R., “Systematic of Vibrational Relaxation,” Journal of 
Chemical Physics, Vol. 139, 1963, pp. 3209-3213. 


4 




5 












Figure 3 . Modified Air Spectrum of the Free-Slream Arc-Jet Flow lor Run 4. 
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Figure 4. Nitric Oxide Potentials Curves Shown on an Energy Diagram. 
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